pacman::p_load(sf, raster, spatstat, sparr, tmap, tidyverse, magick, dplyr)Take-home Exercise 1 - Part 1
Harnessing Geospatial Analytics to Uncover Armed Conflict Patterns in Myanmar
1. Overview
1.1 Background of Myanmar’s Long-Standing Conflicts
.png){kind=tracking}
The conflict in Myanmar is not just a result of the coup but is deeply rooted in the country’s decades-old complex ethnic and political landscape, characterised by tensions between the central government and various ethnic minority groups, each with its own armed forces. The post-coup violence has exacerbated these long-standing conflicts, leading to a severe humanitarian crisis, with thousands killed, hundreds of thousands displaced, and widespread human rights abuses reported.
1.2 Objectives of Take-home Exercise 1
As such, Geospatial analytics has become a valuable tool for evaluating and comprehending the intricacies of increasing conflicts. This exercise aims to reveal the spatial and spatio-temporal distribution of armed conflict in Myanmar by leveraging spatial point pattern analysis. Additionally, it aims to gain clearer insights into the geographical and logistical patterns of violence throughout the nation.
By the end of this take-home exercise, I aim to complete these steps in my spatial point pattern analysis in uncovering the distribution of armed conflict in Myanmar.
- Using appropriate function of sf and tidyverse packages, import and transform the downloaded armed conflict data and administrative boundary data into sf tibble data.frames.
- Using the geospatial data sets prepared, derive quarterly KDE layers.
- Using the geospatial data sets prepared, perform 2nd-Order Spatial Point Patterns Analysis.
- Using the geospatial data sets prepared, derive quarterly spatio-temporal KDE layers.
- Using the geospatial data sets prepared, perform 2nd-Order Spatio-temporal Point Patterns Analysis.
- Using appropriate tmap functions, display the KDE and Spatio-temporal KDE layers on openstreetmap of Myanmar.
- Describe the spatial patterns revealed by the KDE and Spatio-temporal KDE maps.
1.3 About the Datasets
1) Armed Conflict Data (From ACLED)
This Armed Conflict Location & Event Data (ACLED) is an independent, impartial, international non-profit organisation which owns an extensive database of violent conflict and protest in countries and territories around the world.
For the purpose of this exercise, I have downloaded ACLED’s data on Myanmar which includes a series of conflict events, particularly between 1 January 2021 to 30 June 2024.
🔗 Source: ACLED
📁 Format: comma separated values (CSV)
As the dataset is rather extensive, I will be performing my analysis on armed conflict events in a quarterly basis to streamline my tasks. The data included in this dataset are as follows:
Event Type ACLED categorises events into various types. I will mainly be focusing on these four event types: Battles, Explosion/Remote violence, Strategic developments, and Violence against civilians.
|
Location and Geospatial Data The database provides detailed geographic information, pinpointing the exact or approximate locations of conflict events across Myanmar. This includes cities, towns, and rural areas.
|
Date and Time ACLED records the specific dates and, where possible, times of conflict events.
|
Actors
|
Fatalities
|
Others
|
2) Geospatial Data (From Myanmar Information Management Unit)
I will also be using a geospatial dataset from the Myanmar Information Management Unit (MIMU) in shapefile (.shp) format, specifically of the Myanmar state at the 2nd administrative level with district boundaries.
🔗 Source: MIMU
📁 Format: shapefile (.shp)
My reason for choosing the district boundary dataset is that we do not want to select a boundary dataset that is too generalised when analysing conflict events since it might not provide sufficient insights to trends where conflict events happen. Neither do we want to analyse a geography that is too detailed (e.g. Admin 3 - townships) since it can be computationally inefficient as seen in the types of boundary data below.
I have donwloaded the two data sets and organised them into my folder as follows.
2. Let’s Set Up!
2.1 Importing Libraries into R
To carry out this exercise, I will be using the following R packages:
- sf: a relatively new R package specially designed to import, manage and process vector-based geospatial data in R.
- spatstat: has a wide range of useful functions for point pattern analysis. In this take-home exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.
- raster: reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this take-home exercise, it will be used to convert image output generate by spatstat into raster format.
- tmap / sparr: provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.
- tidyverse / dplyr: for transforming and organising data for analysis
- magick: used for plotting animated map plots into GIFs for the spatio-temporal point pattern analysis
Now, let’s install and load these packages in RStudio.
2.2 Importing Data Sets into R
1) Armed Conflicts Data
Next, I will import the downloaded armed conflict data. For aspatial datasets like this, we will import into Rstudio using read_csv() function of the readr package.
# Import armed conflict data
conflict_data <- read_csv("data/aspatial/2021-01-01-2024-06-30-Myanmar.csv")Rows: 87746 Columns: 28
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (18): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (10): year, time_precision, inter1, interaction, iso, latitude, longitud...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Observations
The 2021-01-01-2024-06-30-Myanmar.csv dataset contains 87746 rows and 28 columns which indicates the presence of 87746 unique armed conflict events in Myanmar.
After importing the dataset, we can inspect the dataset using the glimpse() function.
# Inspect the conflict data
glimpse(conflict_data)Rows: 87,746
Columns: 28
$ event_id_cnty <chr> "MMR64313", "MMR64313", "MMR64320", "MMR64320", "MM…
$ event_date <chr> "30 June 2024", "30 June 2024", "30 June 2024", "30…
$ year <dbl> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 202…
$ time_precision <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ disorder_type <chr> "Political violence", "Political violence", "Politi…
$ event_type <chr> "Battles", "Battles", "Battles", "Battles", "Battle…
$ sub_event_type <chr> "Armed clash", "Armed clash", "Armed clash", "Armed…
$ actor1 <chr> "People's Defense Force - Mandalay", "Military Forc…
$ assoc_actor_1 <chr> "MDA - AGF: Madaya - The Authentic Genes Force; SST…
$ inter1 <dbl> 3, 1, 3, 1, 3, 1, 1, 3, 1, 1, 1, 2, 2, 1, 1, 2, 1, …
$ interaction <dbl> 13, 13, 13, 13, 13, 13, 10, 13, 13, 10, 12, 12, 12,…
$ civilian_targeting <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ iso <dbl> 104, 104, 104, 104, 104, 104, 104, 104, 104, 104, 1…
$ region <chr> "Southeast Asia", "Southeast Asia", "Southeast Asia…
$ country <chr> "Myanmar", "Myanmar", "Myanmar", "Myanmar", "Myanma…
$ admin1 <chr> "Mandalay", "Mandalay", "Mandalay", "Mandalay", "Ma…
$ admin2 <chr> "Mandalay", "Mandalay", "Pyinoolwin", "Pyinoolwin",…
$ admin3 <chr> "Patheingyi", "Patheingyi", "Singu", "Singu", "Thab…
$ location <chr> "Aung Tha Pyay", "Aung Tha Pyay", "Pin Lel Gyi", "P…
$ latitude <dbl> 22.1504, 22.1504, 22.5752, 22.5752, 22.8800, 22.880…
$ longitude <dbl> 96.2364, 96.2364, 96.0661, 96.0661, 95.9700, 95.970…
$ geo_precision <dbl> 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 1, …
$ source <chr> "Democratic Voice of Burma; Irrawaddy", "Democratic…
$ source_scale <chr> "National", "National", "National", "National", "Na…
$ notes <chr> "On 30 June 2024, near Aung Tha Pyay village (Pathe…
$ fatalities <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 0, …
$ tags <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ timestamp <dbl> 1720552468, 1720552468, 1720552468, 1720552468, 172…
Observations
The event_date field shows that it uses a character datatype instead of date - we will fix this later. Also, we can observe that thelongitude and langitude fields appear to be adopting the WGS84 geographic coordinate system since they are in the -180/180 and -90/90 range respectively.
2) Myanmar Boundary Data
Observations
When working with Myanmar’s boundary, we need to assign the appropriate coordinate reference system. However, since Myanmar is split into two UTM - West Myanmar (crs: 32646) and East Myanmar (crs: 32647).
Hence, I will also import the administrative boundary data into a simple features tibble data.frame using st_read() of the sf package and check the number of rows returned for both CRS 32646 and 32647. This function reads the shapefile data and returns an sf object that can be used for further analysis.
Find out conflicts count by CRS
conflict_crs <- st_as_sf(conflict_data, coords = c("longitude", "latitude"), crs = 4326)
# Count number of conflicts for CRS 32646
conflict_data_32646 <- st_transform(conflict_crs, crs = 32646)
count_32646 <- nrow(conflict_data_32646)
# Count number of conflicts for CRS 32647
conflict_data_32647 <- st_transform(conflict_crs, crs = 32647)
count_32647 <- nrow(conflict_data_32647)
crs_counts <- data.frame(
CRS = c("EPSG: 32646", "EPSG: 32647"),
Conflicts_Count = c(count_32646, count_32647)
)
print(crs_counts) CRS Conflicts_Count
1 EPSG: 32646 87746
2 EPSG: 32647 87746Since there is no difference in the count, I will decide to focus on UTM zone 47N (EPSG:32647), east of Myanmar, for the purpose of this exercise. The st_transform() function below converts the CRS of the sf object to EPSG:32647.
# Import boundary data
boundary_sf <- st_read(dsn = "data/geospatial",layer = "mmr_polbnda_adm2_250k_mimu") %>% st_transform(crs = 32647)Reading layer `mmr_polbnda_adm2_250k_mimu' from data source
`C:\SamanthaxFoo\IS415-GAA\Take-home_Ex\Take-home_Ex1\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84In the code below, we can notice that the ESPG code has been updated to 32647.
# Check for changes
st_crs(boundary_sf)Coordinate Reference System:
User input: EPSG:32647
wkt:
PROJCRS["WGS 84 / UTM zone 47N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 47N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",99,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
BBOX[0,96,84,102]],
ID["EPSG",32647]]Here, I will use the plot() function which plots the geometry of the sf object. The st_geometry() function is used to extract the geometry of the mpsz_sf object which includes the districts of Myanmar as shown below.
par(mar = c(0,0,0,0))
plot(st_geometry(boundary_sf))
3. Data Wrangling
3.1 Fixing Incorrect Datatypes
Recall that the earlier inspection of the conflict_data tibble data frame revealed that the datatype indicated for event date is wrongly labelled as a character instead of a date format.
As such, let’s convert the datatype to the correct ‘date’ format as shown below.
# Convert the datatype for event_date
conflict_data$event_date <- as.Date(conflict_data$event_date, format = "%d %B %Y")
# Check for changes
head(conflict_data)# A tibble: 6 × 28
event_id_cnty event_date year time_precision disorder_type event_type
<chr> <date> <dbl> <dbl> <chr> <chr>
1 MMR64313 2024-06-30 2024 1 Political violence Battles
2 MMR64313 2024-06-30 2024 1 Political violence Battles
3 MMR64320 2024-06-30 2024 1 Political violence Battles
4 MMR64320 2024-06-30 2024 1 Political violence Battles
5 MMR64321 2024-06-30 2024 1 Political violence Battles
6 MMR64321 2024-06-30 2024 1 Political violence Battles
# ℹ 22 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
# inter1 <dbl>, interaction <dbl>, civilian_targeting <chr>, iso <dbl>,
# region <chr>, country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>,
# location <chr>, latitude <dbl>, longitude <dbl>, geo_precision <dbl>,
# source <chr>, source_scale <chr>, notes <chr>, fatalities <dbl>,
# tags <chr>, timestamp <dbl>3.2 Adding new year_quarter column
We will want to create a new column to indicate the specific year and quarter for each conflict event since the spatial analysis will be done later in a quarterly manner.
Extract year and quarter
conflict_data$year_quarter <- paste0(
year(conflict_data$event_date),
" Q",
quarter(conflict_data$event_date)
)
# View the new data column
unique(conflict_data$year_quarter) [1] "2024 Q2" "2024 Q1" "2023 Q4" "2023 Q3" "2023 Q2" "2023 Q1" "2022 Q4"
[8] "2022 Q3" "2022 Q2" "2022 Q1" "2021 Q4" "2021 Q3" "2021 Q2" "2021 Q1"3.3 Fixing Duplicated Event ID in conflict_data Dataframe
As shown, there are presence of duplicates in our dataframe returned by the duplicated() function.
# Check for duplicates
any(duplicated(conflict_data))[1] TRUEBased on the duplicated event ID: MMR64313 for instance. We can observe the two records are of the same political violence event happening between two actors on 30/6/2024, between the People’s Defense Force and Military Forces of Myanmar. Upon further research, these two actors are opposing political parties of Myanmar’s ongoing conflict.
# Inspect an instance of the duplciated event IDs
head(conflict_data,2)# A tibble: 2 × 29
event_id_cnty event_date year time_precision disorder_type event_type
<chr> <date> <dbl> <dbl> <chr> <chr>
1 MMR64313 2024-06-30 2024 1 Political violence Battles
2 MMR64313 2024-06-30 2024 1 Political violence Battles
# ℹ 23 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
# inter1 <dbl>, interaction <dbl>, civilian_targeting <chr>, iso <dbl>,
# region <chr>, country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>,
# location <chr>, latitude <dbl>, longitude <dbl>, geo_precision <dbl>,
# source <chr>, source_scale <chr>, notes <chr>, fatalities <dbl>,
# tags <chr>, timestamp <dbl>, year_quarter <chr>
Reflection
Should duplicated data be removed in this analysis?
A single event (e.g. MMR64313) can have duplicated rows with different actor1 values, typically due to counterattacks from opposing sides, leading to different data entries into the conflict_data dataset.
Hence, I will remove duplicated events found in the conflict_data dataframe as long as the rows have the same event ID indicated.
Here, I did another check to ensure there is not more than 2 possible repeated event IDs in the first 20 rows of conflict_data.
Check duplicated events for first 20 rows
duplicate_counts_first_20 <- conflict_data %>%
slice(1:20) %>%
group_by(event_id_cnty) %>%
summarize(count = n()) %>%
filter(count > 1)
# View the result
print(duplicate_counts_first_20)# A tibble: 9 × 2
event_id_cnty count
<chr> <int>
1 MMR64313 2
2 MMR64320 2
3 MMR64321 2
4 MMR64323 2
5 MMR64325 2
6 MMR64326 2
7 MMR64328 2
8 MMR64330 2
9 MMR64331 2With that checked, I’ll remove the duplicated rows with a repeated Event ID but will retain the actor1 value and append it to a new actor2 column, assuming two actors are involved in every conflict.
Remove duplicated rows
# Load necessary library
library(dplyr)
# Retrieve data of duplicated rows and summarize actor2, assoc_actor_2, and inter2
merged_duplicates <- conflict_data %>%
filter(duplicated(event_id_cnty) | duplicated(event_id_cnty, fromLast = TRUE)) %>%
arrange(event_id_cnty) %>%
group_by(event_id_cnty) %>%
summarize(
actor2 = last(actor1)
)
# Keep rows without duplicates
conflict_data_no_duplicates <- conflict_data %>%
filter(!duplicated(event_id_cnty))
# Update conflict_data dataframe with new columns from merged_duplicates
conflict_data <- conflict_data_no_duplicates %>%
left_join(merged_duplicates, by = "event_id_cnty") %>%
mutate(
# Ensure actor2 exists after the join, if not fill with actor1
actor2 = coalesce(actor2, actor1)
)
# View dataframe
head(conflict_data[c('event_id_cnty','actor1','actor2')])# A tibble: 6 × 3
event_id_cnty actor1 actor2
<chr> <chr> <chr>
1 MMR64313 People's Defense Force - Mandalay Military Force…
2 MMR64320 People's Defense Force - Mandalay Military Force…
3 MMR64321 People's Defense Force - Mandalay Military Force…
4 MMR64322 Military Forces of Myanmar (2021-) Military Force…
5 MMR64323 PKDF MMU: People's Knight Defense Force - Myinmu Military Force…
6 MMR64324 Military Forces of Myanmar (2021-) Military Force…We can observe that there are no longer any duplicated event IDs in our conflict_data data frame.
any(duplicated(conflict_data))[1] FALSE3.4 Converting Aspatial Data to Simple Feature Format
For the purpose of this exercise, we will want to integrate and analyse aspatial data in a geographic context. I’ll do a check if conflict_data needs to be converted to a sf data frame - if it outputs anything else but sf, then it’s not a simple feature data frame!
class(conflict_data)[1] "tbl_df" "tbl" "data.frame"
Observations
We can see that conflict_data is not a sf data frame. Since a non-simple feature data frame does not have a “geometry” column, we’ll need to convert conflict_data into a simple feature data frame
We can convert conflict_data into a simple feature data frame by using st_as_sf() from the sf package. Addiitionally, we will also need to transform coordinate system from geographic (ESPG: 4326) to projected (ESPG: 32647) using st_transform().
# Convert to simple feature format
conflict_data_sf <- st_as_sf(conflict_data, coords = c("longitude", "latitude"), crs=4326) %>% st_transform(crs = 32647)
# Inspect the changes
glimpse(conflict_data_sf)Rows: 51,553
Columns: 29
$ event_id_cnty <chr> "MMR64313", "MMR64320", "MMR64321", "MMR64322", "MM…
$ event_date <date> 2024-06-30, 2024-06-30, 2024-06-30, 2024-06-30, 20…
$ year <dbl> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 202…
$ time_precision <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ disorder_type <chr> "Political violence", "Political violence", "Politi…
$ event_type <chr> "Battles", "Battles", "Battles", "Strategic develop…
$ sub_event_type <chr> "Armed clash", "Armed clash", "Armed clash", "Chang…
$ actor1 <chr> "People's Defense Force - Mandalay", "People's Defe…
$ assoc_actor_1 <chr> "MDA - AGF: Madaya - The Authentic Genes Force; SST…
$ inter1 <dbl> 3, 3, 3, 1, 3, 1, 1, 2, 1, 1, 1, 1, 3, 3, 3, 7, 1, …
$ interaction <dbl> 13, 13, 13, 10, 13, 10, 12, 12, 12, 12, 12, 13, 13,…
$ civilian_targeting <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ iso <dbl> 104, 104, 104, 104, 104, 104, 104, 104, 104, 104, 1…
$ region <chr> "Southeast Asia", "Southeast Asia", "Southeast Asia…
$ country <chr> "Myanmar", "Myanmar", "Myanmar", "Myanmar", "Myanma…
$ admin1 <chr> "Mandalay", "Mandalay", "Mandalay", "Sagaing", "Sag…
$ admin2 <chr> "Mandalay", "Pyinoolwin", "Pyinoolwin", "Shwebo", "…
$ admin3 <chr> "Patheingyi", "Singu", "Thabeikkyin", "Khin-U", "My…
$ location <chr> "Aung Tha Pyay", "Pin Lel Gyi", "Thabeikkyin", "Khi…
$ geo_precision <dbl> 2, 1, 2, 2, 2, 2, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 1, …
$ source <chr> "Democratic Voice of Burma; Irrawaddy", "Irrawaddy"…
$ source_scale <chr> "National", "National", "National", "Subnational-Na…
$ notes <chr> "On 30 June 2024, near Aung Tha Pyay village (Pathe…
$ fatalities <dbl> 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ tags <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ timestamp <dbl> 1720552468, 1720552468, 1720552468, 1720552468, 172…
$ year_quarter <chr> "2024 Q2", "2024 Q2", "2024 Q2", "2024 Q2", "2024 Q…
$ actor2 <chr> "Military Forces of Myanmar (2021-)", "Military For…
$ geometry <POINT [m]> POINT (214961 2452068), POINT (198303.2 24994…
Observations
Notice that a new column called geometry has been added into the data frame. On the other hand, the longitude and latitude columns have been removed from the data frame.
We can further inspect the newly created ‘geometry’ column of conflict_data_sf
# Retrieve geometry column
st_geometry(conflict_data_sf)Geometry set for 51553 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -208804.4 ymin: 1103500 xmax: 640934.5 ymax: 3042960
Projected CRS: WGS 84 / UTM zone 47N
First 5 geometries:POINT (214961 2452068)POINT (198303.2 2499463)POINT (189105.4 2533434)POINT (160913.9 2522331)POINT (146213 2428487)
Observations
It consists of 51,533 features consisting of point geometric features where the underlying datum is in WGS 84 format.
To ensure that the coordinate system is correctly updated, we can use the st_crs() function where we observe that the ESPG code is correctly indicated as 32647.
# Check CRS format
st_crs(conflict_data_sf)Coordinate Reference System:
User input: EPSG:32647
wkt:
PROJCRS["WGS 84 / UTM zone 47N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 47N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",99,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
BBOX[0,96,84,102]],
ID["EPSG",32647]]3.5 Reduce Data File Size
In this section, I will reduce the current Myanmar armed conflict dataset as the time taken for computing the kernel density estimates can take up to 30 minutes long which is not computationally efficient.
1) Remove ‘Protests’ and ‘Riots’ Event Types
I will remove rows in the conflicts_data_sf dataset that don’t focus on the four main event types (Battles, Explosion/Remote violence, Strategic developments, and Violence against civilians), as mentioned in the exercise brief.
conflict_data_sf <- conflict_data_sf %>%
filter(!(event_type %in% c("Protests", "Riots")))
unique(conflict_data_sf$event_type)[1] "Battles" "Strategic developments"
[3] "Violence against civilians" "Explosions/Remote violence"2) Remove unused columns in boundary_sf
As seen, there are 8 columns in the simple feature data frame of boundary_sf.
# Inspect first rows of data in boundary_sf
head(boundary_sf)Simple feature collection with 6 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -14915.04 ymin: 1736124 xmax: 187961.7 ymax: 2051144
Projected CRS: WGS 84 / UTM zone 47N
OBJECTID ST ST_PCODE DT DT_PCODE DT_MMR PCode_V
1 1 Ayeyarwady MMR017 Hinthada MMR017D002 ဟင်္သာတခရိုင် 9.4
2 2 Ayeyarwady MMR017 Labutta MMR017D004 လပွတ္တာခရိုင် 9.4
3 3 Ayeyarwady MMR017 Maubin MMR017D005 မအူပင်ခရိုင် 9.4
4 4 Ayeyarwady MMR017 Myaungmya MMR017D003 မြောင်းမြခရိုင် 9.4
5 5 Ayeyarwady MMR017 Pathein MMR017D001 ပုသိမ်ခရိုင် 9.4
6 6 Ayeyarwady MMR017 Pyapon MMR017D006 ဖျာပုံခရိုင် 9.4
geometry
1 MULTIPOLYGON (((90859.89 20...
2 MULTIPOLYGON (((75991.51 17...
3 MULTIPOLYGON (((115559 1928...
4 MULTIPOLYGON (((39919.39 18...
5 MULTIPOLYGON (((-6302.348 1...
6 MULTIPOLYGON (((93411.72 17...I will remove ’DT_MMR” column as we already have the District Name in English in DT and won’t require the district names in Myanmar Language. Next, we will remove the coded versions of ST (state/region) and DT (district) columns, namely ST_PCODE and DT_PCODE. Additionally, we won’t need the PCode_V column since we will be dropping the PCODE column too.
boundary_sf <- boundary_sf %>% dplyr::select('OBJECTID', 'ST', 'DT','geometry')
summary(boundary_sf) OBJECTID ST DT geometry
Min. : 1.00 Length:80 Length:80 MULTIPOLYGON :80
1st Qu.:20.75 Class :character Class :character epsg:32647 : 0
Median :40.50 Mode :character Mode :character +proj=utm ...: 0
Mean :40.50
3rd Qu.:60.25
Max. :80.00 3) Remove unused columns in conflict_data
I will also remove unnecessary columns of the conflict_data data frame that won’t be used in our spatial analysis later. I’ll rename admin1 to ST (states) and admin2 to DT (districts) for easier reference.
Remove unnecessary columns
conflict_data_sf <- conflict_data_sf %>%
select(event_id_cnty, event_date, year_quarter, disorder_type,
event_type, location, admin1, admin2, geometry, actor1, actor2,
interaction, fatalities) %>%
rename(
ST = admin1,
DT = admin2
)
summary(conflict_data_sf) event_id_cnty event_date year_quarter disorder_type
Length:42608 Min. :2021-01-01 Length:42608 Length:42608
Class :character 1st Qu.:2022-01-10 Class :character Class :character
Mode :character Median :2022-10-13 Mode :character Mode :character
Mean :2022-10-29
3rd Qu.:2023-08-29
Max. :2024-06-30
event_type location ST DT
Length:42608 Length:42608 Length:42608 Length:42608
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
geometry actor1 actor2 interaction
POINT :42608 Length:42608 Length:42608 Min. :10.00
epsg:32647 : 0 Class :character Class :character 1st Qu.:13.00
+proj=utm ...: 0 Mode :character Mode :character Median :17.00
Mean :18.86
3rd Qu.:17.00
Max. :80.00
fatalities
Min. : 0.00
1st Qu.: 0.00
Median : 0.00
Mean : 1.27
3rd Qu.: 1.00
Max. :201.00 3.6 Converting Simple Features Data Frame into ppp Object
It is important that we convert conflict_data_sf (a simple feature data frame) into a planer point pattern (ppp) object format, since the spatstat package that we’ll be using for the Spatial Point Pattern Analysis later is specifically designed for working with ppp-formated data. Additionally, I will begin with categorising the ppp objects into their unique year_quarter category.
Create ppp objects based on year_quarter category
# Create an empty list to store the ppp objects
ppp_list <- list()
# Loop through each unique year_quarter category
for (yq in unique(conflict_data_sf$year_quarter)) {
# Subset the data for the current year_quarter
subset_data_sf <- conflict_data_sf %>% filter(year_quarter == yq)
# Convert the subset to a ppp object
subset_ppp <- as.ppp(subset_data_sf$geometry)
# Add the ppp object to the list
ppp_list[[yq]] <- subset_ppp
}
# Check list
ppp_list$`2024 Q2`
Planar point pattern: 2788 points
window: rectangle = [-208804.4, 597543.7] x [1103500.1, 3026504.9] units
$`2024 Q1`
Planar point pattern: 3186 points
window: rectangle = [-207135, 591875.9] x [1245380, 3026504.9] units
$`2023 Q4`
Planar point pattern: 3627 points
window: rectangle = [-206931.7, 604775.1] x [1103500.1, 3020772.2] units
$`2023 Q3`
Planar point pattern: 3010 points
window: rectangle = [-197883.4, 518300.4] x [1103500.1, 3027041.8] units
$`2023 Q2`
Planar point pattern: 2745 points
window: rectangle = [-191261.5, 518300.4] x [1103500.1, 3006372.9] units
$`2023 Q1`
Planar point pattern: 3101 points
window: rectangle = [-199243.8, 591875.9] x [1103500.1, 3026504.9] units
$`2022 Q4`
Planar point pattern: 3296 points
window: rectangle = [-206531.5, 518300.4] x [1103500.1, 2931517.1] units
$`2022 Q3`
Planar point pattern: 3486 points
window: rectangle = [-206196.6, 568361.5] x [1103500.1, 3026504.9] units
$`2022 Q2`
Planar point pattern: 3580 points
window: rectangle = [-206931.7, 640934.5] x [1103500.1, 3026504.9] units
$`2022 Q1`
Planar point pattern: 3563 points
window: rectangle = [-204784, 591875.9] x [1103500.1, 3026504.9] units
$`2021 Q4`
Planar point pattern: 3844 points
window: rectangle = [-200024.3, 591875.9] x [1103500.1, 3042960.3] units
$`2021 Q3`
Planar point pattern: 2754 points
window: rectangle = [-193181.1, 591875.9] x [1103500.1, 3042960.3] units
$`2021 Q2`
Planar point pattern: 2916 points
window: rectangle = [-191409.1, 640934.5] x [1132472.1, 3042960.3] units
$`2021 Q1`
Planar point pattern: 712 points
window: rectangle = [-203795.3, 591875.9] x [1375186.1, 3026504.9] unitsWe can visualise the spread of conflict events across each quarter from January 2021 to June 2024 using the plot() function as shown below.
Visualise the spread of conflicts by year_quarter
# Ensure 'year_quarter' is a factor
conflict_data_sf$year_quarter <- as.factor(conflict_data_sf$year_quarter)
# Loop through each unique year_quarter and create separate plots
year_quarters <- unique(conflict_data_sf$year_quarter)
# Set up a grid layout for multiple plots (adjust 'mfrow' as needed)
par(mfrow = c(4,4))
par(mar = c(0,0,1,0))
# Loop through each year_quarter and plot
for (yq in year_quarters) {
subset_data_sf <- conflict_data_sf[conflict_data_sf$year_quarter == yq, ]
conflict_data_ppp <- as.ppp(subset_data_sf$geometry)
# Plot each subset ppp object
plot(conflict_data_ppp, main = paste(yq))
}
Observations
It is noticeable that there conflict events have occured more frequently since 2021 as points plotted on the graph have gotten darker across 2021 to 2024. We can also observe the possibility of duplicated events occurring from the darker spots in the plot, in which it appears more intense in Myanmar’s central and west regions.
3.7 Creating owin object
When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area, that is Myanmar’s boundary in this case. In spatstat, an object called owin is specially designed to represent this polygonal region.
The code chunk below is used to convert the boundary_data_sf simple feature data frame into an owin object of spatstat.
# Convert to owin object
myanmar_owin <- as.owin(boundary_sf)
# Visualise the owin object
plot(myanmar_owin)
Observations
From my observations, the as.owin() function converts the boundary_data_sf spatial boundary into a window object that represents the outer boundary of the spatial region and does not handle internal structures or districts we previously saw from the plot of boundary_data_sf.
We can also take a quick look at the owin object properties as shown. I will be converting it to a data frame for the purposes of getting a quick glimpse of the object.
# Summary info of owin object
owin_df <- as.data.frame(myanmar_owin)
print(head(owin_df)) x y id sign
1 56519.39 2741919 1 -1
2 56917.28 2741947 1 -1
3 57000.15 2741973 1 -1
4 57068.51 2741994 1 -1
5 57221.44 2742142 1 -1
6 57068.51 2741994 1 -13.8 Combining ppp Object and owin Object
In this last step of geospatial data wrangling, I will mask all ppp object with the owin object I created earlier to put in place all conflict events within the boundary of Myanmar. Doing so can also optimise the memory usage for large datasets. I’ll also make the density values more comprehensible by re-scaling the density values from metres to kilometres using rescale().
# Initialize an empty list to store masked ppp objects
masked_ppp_list_km <- list()
# Iterate over each ppp object in the list
for (quarter in names(ppp_list)) {
ppp_obj <- ppp_list[[quarter]]
# Mask the ppp object with the owin object
masked_ppp <- ppp_obj[myanmar_owin]
masked_ppp_km <- rescale(masked_ppp, 1000, "km")
# Store the masked ppp object in the new list
masked_ppp_list_km[[quarter]] <- masked_ppp_km
}# Inspect 2024 Q2 masked ppp object
summary(masked_ppp_list_km$`2024 Q2`)Planar point pattern: 2788 points
Average intensity 0.004162974 points per square km
*Pattern contains duplicated points*
Coordinates are given to 16 decimal places
Window: polygonal boundary
1077 separate polygons (510 holes)
vertices area relative.area
polygon 1 (hole) 9 -4.40155e-08 -6.57e-14
polygon 2 (hole) 3 -4.82818e-08 -7.21e-14
polygon 3 (hole) 3 -2.78345e-09 -4.16e-15
polygon 4 (hole) 3 -4.99415e-09 -7.46e-15
polygon 5 (hole) 3 -6.70916e-10 -1.00e-15
polygon 6 (hole) 3 -4.99758e-09 -7.46e-15
polygon 7 (hole) 6 -3.24429e-09 -4.84e-15
polygon 8 (hole) 11 -5.35507e-08 -8.00e-14
polygon 9 (hole) 3 -3.06875e-09 -4.58e-15
polygon 10 (hole) 3 -7.63672e-09 -1.14e-14
polygon 11 (hole) 3 -8.61828e-11 -1.29e-16
polygon 12 (hole) 3 -2.44104e-08 -3.64e-14
polygon 13 (hole) 3 -4.88141e-11 -7.29e-17
polygon 14 (hole) 4 -1.74466e-10 -2.61e-16
polygon 15 (hole) 4 -2.30890e-11 -3.45e-17
polygon 16 (hole) 4 -6.25346e-14 -9.34e-20
polygon 17 (hole) 4 -3.68040e-08 -5.50e-14
polygon 18 (hole) 3 -1.27228e-12 -1.90e-18
polygon 19 (hole) 4 -1.25504e-09 -1.87e-15
polygon 20 (hole) 3 -1.56980e-08 -2.34e-14
polygon 21 (hole) 4 -1.13609e-12 -1.70e-18
polygon 22 (hole) 4 -8.02868e-08 -1.20e-13
polygon 23 (hole) 4 -3.98057e-11 -5.94e-17
polygon 24 (hole) 4 -2.89739e-08 -4.33e-14
polygon 25 (hole) 5 -7.39370e-08 -1.10e-13
polygon 26 (hole) 6 -1.08246e-07 -1.62e-13
polygon 27 (hole) 13 -3.27066e-07 -4.88e-13
polygon 28 (hole) 3 -3.55050e-09 -5.30e-15
polygon 29 (hole) 3 -2.84064e-08 -4.24e-14
polygon 30 (hole) 3 -8.07641e-10 -1.21e-15
polygon 31 (hole) 4 -5.20948e-09 -7.78e-15
polygon 32 (hole) 4 -3.21399e-10 -4.80e-16
polygon 33 (hole) 4 -2.54955e-11 -3.81e-17
polygon 34 (hole) 3 -4.29638e-10 -6.42e-16
polygon 35 (hole) 6 -6.02383e-10 -8.99e-16
polygon 36 (hole) 4 -4.91852e-10 -7.34e-16
polygon 37 26 2.85778e+00 4.27e-06
polygon 38 (hole) 3 -6.80850e-13 -1.02e-18
polygon 39 (hole) 4 -3.09085e-13 -4.62e-19
polygon 40 (hole) 6 -3.61049e-12 -5.39e-18
polygon 41 (hole) 7 -5.33310e-12 -7.96e-18
polygon 42 (hole) 4 -3.62663e-12 -5.42e-18
polygon 43 (hole) 3 -2.00302e-14 -2.99e-20
polygon 44 (hole) 3 -2.85255e-12 -4.26e-18
polygon 45 (hole) 4 -1.15535e-12 -1.73e-18
polygon 46 (hole) 17 -6.25657e-13 -9.34e-19
polygon 47 (hole) 20 -1.53263e-11 -2.29e-17
polygon 48 (hole) 3 -1.42406e-14 -2.13e-20
polygon 49 3 0.00000e+00 0.00e+00
polygon 50 (hole) 4 -2.46048e-12 -3.67e-18
polygon 51 (hole) 4 -2.93672e-13 -4.39e-19
polygon 52 (hole) 3 -1.35167e-16 -2.02e-22
polygon 53 (hole) 3 -3.66243e-13 -5.47e-19
polygon 54 (hole) 13 -3.22774e-12 -4.82e-18
polygon 55 (hole) 7 -3.79713e-12 -5.67e-18
polygon 56 (hole) 4 -6.06083e-13 -9.05e-19
polygon 57 (hole) 6 -5.94835e-12 -8.88e-18
polygon 58 (hole) 4 -4.11305e-13 -6.14e-19
polygon 59 (hole) 12 -7.16981e-12 -1.07e-17
polygon 60 (hole) 4 -1.96241e-14 -2.93e-20
polygon 61 (hole) 4 -6.16304e-12 -9.20e-18
polygon 62 (hole) 4 -2.13366e-13 -3.19e-19
polygon 63 (hole) 4 -2.21852e-12 -3.31e-18
polygon 64 (hole) 3 -4.22166e-15 -6.30e-21
polygon 65 (hole) 5 -6.21282e-13 -9.28e-19
polygon 66 (hole) 4 -5.02821e-16 -7.51e-22
polygon 67 (hole) 4 -3.30895e-13 -4.94e-19
polygon 68 (hole) 17 -1.19462e-11 -1.78e-17
polygon 69 (hole) 3 -1.27084e-13 -1.90e-19
polygon 70 (hole) 4 -1.80406e-13 -2.69e-19
polygon 71 (hole) 8 -3.87178e-12 -5.78e-18
polygon 72 (hole) 4 -7.14736e-13 -1.07e-18
polygon 73 (hole) 5 -3.43923e-12 -5.14e-18
polygon 74 3 0.00000e+00 0.00e+00
polygon 75 (hole) 8 -2.57732e-12 -3.85e-18
polygon 76 (hole) 4 -1.36629e-12 -2.04e-18
polygon 77 (hole) 13 -1.46870e-11 -2.19e-17
polygon 78 (hole) 3 -3.08778e-14 -4.61e-20
polygon 79 (hole) 18 -1.78311e-12 -2.66e-18
polygon 80 (hole) 28 -1.67874e-11 -2.51e-17
polygon 81 (hole) 4 -4.46113e-13 -6.66e-19
polygon 82 (hole) 3 -8.98160e-15 -1.34e-20
polygon 83 (hole) 4 -1.43028e-12 -2.14e-18
polygon 84 (hole) 4 -6.38438e-14 -9.53e-20
polygon 85 (hole) 19 -1.89999e-11 -2.84e-17
polygon 86 (hole) 4 -1.49124e-13 -2.23e-19
polygon 87 (hole) 3 -3.68533e-14 -5.50e-20
polygon 88 (hole) 4 -4.45374e-13 -6.65e-19
polygon 89 (hole) 4 -4.83479e-13 -7.22e-19
polygon 90 (hole) 12 -9.32407e-12 -1.39e-17
polygon 91 (hole) 3 -8.84897e-13 -1.32e-18
polygon 92 (hole) 12 -8.41711e-12 -1.26e-17
polygon 93 (hole) 3 -3.73405e-14 -5.58e-20
polygon 94 (hole) 3 -1.40773e-12 -2.10e-18
polygon 95 (hole) 4 -7.84182e-13 -1.17e-18
polygon 96 (hole) 5 -4.10716e-13 -6.13e-19
polygon 97 (hole) 10 -5.89467e-12 -8.80e-18
polygon 98 (hole) 3 -3.54370e-15 -5.29e-21
polygon 99 (hole) 11 -6.90903e-12 -1.03e-17
polygon 100 (hole) 3 -1.89866e-13 -2.84e-19
polygon 101 (hole) 3 -3.09480e-15 -4.62e-21
polygon 102 (hole) 3 -3.95637e-17 -5.91e-23
polygon 103 (hole) 4 -1.88731e-16 -2.82e-22
polygon 104 (hole) 3 -5.93215e-17 -8.86e-23
polygon 105 (hole) 6 -1.10389e-12 -1.65e-18
polygon 106 (hole) 5 -1.62107e-13 -2.42e-19
polygon 107 (hole) 4 -6.07269e-13 -9.07e-19
polygon 108 (hole) 3 -4.17062e-14 -6.23e-20
polygon 109 (hole) 7 -4.10998e-12 -6.14e-18
polygon 110 (hole) 7 -2.99640e-12 -4.47e-18
polygon 111 (hole) 4 -8.08268e-13 -1.21e-18
polygon 112 (hole) 8 -1.72022e-12 -2.57e-18
polygon 113 (hole) 3 -5.46853e-13 -8.17e-19
polygon 114 (hole) 3 -9.21233e-15 -1.38e-20
polygon 115 (hole) 3 -2.34233e-14 -3.50e-20
polygon 116 (hole) 9 -1.27183e-12 -1.90e-18
polygon 117 (hole) 4 -8.58065e-13 -1.28e-18
polygon 118 (hole) 6 -4.85387e-12 -7.25e-18
polygon 119 (hole) 6 -4.38024e-12 -6.54e-18
polygon 120 (hole) 6 -1.18459e-12 -1.77e-18
polygon 121 (hole) 18 -5.24919e-12 -7.84e-18
polygon 122 3 0.00000e+00 0.00e+00
polygon 123 (hole) 10 -3.29610e-12 -4.92e-18
polygon 124 (hole) 4 -8.31786e-15 -1.24e-20
polygon 125 (hole) 8 -5.92159e-12 -8.84e-18
polygon 126 (hole) 4 -2.92352e-15 -4.37e-21
polygon 127 (hole) 3 -2.22606e-13 -3.32e-19
polygon 128 (hole) 16 -1.60557e-11 -2.40e-17
polygon 129 (hole) 10 -4.89387e-12 -7.31e-18
polygon 130 (hole) 4 -1.48133e-13 -2.21e-19
polygon 131 (hole) 4 -2.04528e-12 -3.05e-18
polygon 132 (hole) 4 -4.07107e-13 -6.08e-19
polygon 133 (hole) 4 -2.13923e-12 -3.19e-18
polygon 134 (hole) 4 -2.58593e-17 -3.86e-23
polygon 135 (hole) 3 -8.40799e-14 -1.26e-19
polygon 136 (hole) 4 -5.03212e-13 -7.51e-19
polygon 137 (hole) 4 -1.10455e-12 -1.65e-18
polygon 138 (hole) 4 -3.74743e-14 -5.60e-20
polygon 139 (hole) 6 -1.42384e-12 -2.13e-18
polygon 140 (hole) 4 -1.96713e-14 -2.94e-20
polygon 141 (hole) 4 -7.29770e-13 -1.09e-18
polygon 142 (hole) 3 -7.99677e-16 -1.19e-21
polygon 143 (hole) 3 -1.00370e-15 -1.50e-21
polygon 144 (hole) 3 -9.07286e-16 -1.35e-21
polygon 145 (hole) 4 -2.63668e-14 -3.94e-20
polygon 146 (hole) 5 -8.51993e-13 -1.27e-18
polygon 147 (hole) 4 -3.15631e-13 -4.71e-19
polygon 148 (hole) 4 -7.89181e-13 -1.18e-18
polygon 149 (hole) 6 -4.24179e-13 -6.33e-19
polygon 150 (hole) 16 -1.00987e-11 -1.51e-17
polygon 151 (hole) 3 -6.53241e-18 -9.75e-24
polygon 152 (hole) 6 -2.87035e-12 -4.29e-18
polygon 153 (hole) 3 -2.64530e-15 -3.95e-21
polygon 154 (hole) 8 -1.99962e-12 -2.99e-18
polygon 155 (hole) 3 -2.45345e-13 -3.66e-19
polygon 156 (hole) 3 -3.99767e-14 -5.97e-20
polygon 157 (hole) 6 -1.08025e-12 -1.61e-18
polygon 158 (hole) 3 -1.22938e-13 -1.84e-19
polygon 159 (hole) 5 -3.07870e-13 -4.60e-19
polygon 160 (hole) 12 -4.81197e-12 -7.19e-18
polygon 161 (hole) 4 -8.16819e-14 -1.22e-19
polygon 162 (hole) 3 -1.28714e-12 -1.92e-18
polygon 163 (hole) 4 -1.87938e-13 -2.81e-19
polygon 164 (hole) 11 -1.79968e-12 -2.69e-18
polygon 165 (hole) 3 -5.93827e-15 -8.87e-21
polygon 166 (hole) 3 -4.50189e-13 -6.72e-19
polygon 167 (hole) 4 -4.87377e-13 -7.28e-19
polygon 168 (hole) 13 -1.01705e-11 -1.52e-17
polygon 169 (hole) 3 -7.37542e-13 -1.10e-18
polygon 170 (hole) 4 -3.62179e-13 -5.41e-19
polygon 171 3 9.13135e-19 1.36e-24
polygon 172 (hole) 8 -1.74082e-12 -2.60e-18
polygon 173 (hole) 3 -9.31514e-14 -1.39e-19
polygon 174 43 7.32477e+00 1.09e-05
polygon 175 (hole) 3 -6.79111e-14 -1.01e-19
polygon 176 (hole) 19 -1.27551e-11 -1.90e-17
polygon 177 (hole) 3 -7.65776e-15 -1.14e-20
polygon 178 (hole) 9 -4.34371e-12 -6.49e-18
polygon 179 (hole) 4 -7.09328e-15 -1.06e-20
polygon 180 (hole) 6 -8.43021e-12 -1.26e-17
polygon 181 (hole) 30 -5.76461e-12 -8.61e-18
polygon 182 (hole) 12 -5.45083e-12 -8.14e-18
polygon 183 (hole) 16 -1.90750e-11 -2.85e-17
polygon 184 (hole) 15 -2.10315e-11 -3.14e-17
polygon 185 (hole) 52 -3.02970e-11 -4.52e-17
polygon 186 103 1.86991e+01 2.79e-05
polygon 187 (hole) 10 -7.40529e-12 -1.11e-17
polygon 188 (hole) 6 -2.03062e-12 -3.03e-18
polygon 189 (hole) 7 -1.02421e-12 -1.53e-18
polygon 190 (hole) 3 -1.32600e-19 -1.98e-25
polygon 191 (hole) 3 -5.60223e-13 -8.37e-19
polygon 192 (hole) 10 -2.37610e-12 -3.55e-18
polygon 193 (hole) 3 -1.90831e-16 -2.85e-22
polygon 194 (hole) 4 -1.76514e-12 -2.64e-18
polygon 195 (hole) 3 -1.37598e-16 -2.05e-22
polygon 196 (hole) 4 -5.28260e-13 -7.89e-19
polygon 197 37 9.32316e+00 1.39e-05
polygon 198 (hole) 4 -2.33117e-13 -3.48e-19
polygon 199 (hole) 4 -6.59209e-13 -9.84e-19
polygon 200 (hole) 4 -5.12582e-13 -7.65e-19
polygon 201 (hole) 4 -4.70065e-13 -7.02e-19
polygon 202 (hole) 3 -2.10303e-13 -3.14e-19
polygon 203 (hole) 6 -1.36386e-12 -2.04e-18
polygon 204 (hole) 4 -2.18244e-13 -3.26e-19
polygon 205 (hole) 4 -3.72475e-13 -5.56e-19
polygon 206 (hole) 4 -8.32022e-13 -1.24e-18
polygon 207 (hole) 4 -1.08153e-12 -1.61e-18
polygon 208 (hole) 3 -5.92333e-17 -8.84e-23
polygon 209 (hole) 7 -6.03242e-12 -9.01e-18
polygon 210 (hole) 6 -2.27756e-12 -3.40e-18
polygon 211 (hole) 4 -1.53368e-12 -2.29e-18
polygon 212 (hole) 3 -2.33223e-14 -3.48e-20
polygon 213 (hole) 6 -8.61964e-12 -1.29e-17
polygon 214 371 2.43869e+02 3.64e-04
polygon 215 (hole) 4 -8.45587e-13 -1.26e-18
polygon 216 297 2.84905e+02 4.25e-04
polygon 217 (hole) 4 -3.40250e-13 -5.08e-19
polygon 218 (hole) 6 -1.39288e-12 -2.08e-18
polygon 219 (hole) 4 -1.27223e-13 -1.90e-19
polygon 220 (hole) 4 -9.59107e-13 -1.43e-18
polygon 221 (hole) 8 -1.39768e-11 -2.09e-17
polygon 222 (hole) 3 -3.34368e-14 -4.99e-20
polygon 223 (hole) 4 -9.58519e-13 -1.43e-18
polygon 224 (hole) 5 -1.32345e-13 -1.98e-19
polygon 225 (hole) 4 -5.98296e-13 -8.93e-19
polygon 226 (hole) 7 -3.05539e-12 -4.56e-18
polygon 227 (hole) 3 -1.80098e-15 -2.69e-21
polygon 228 (hole) 9 -6.38785e-12 -9.54e-18
polygon 229 33 1.68222e+01 2.51e-05
polygon 230 (hole) 8 -1.04272e-11 -1.56e-17
polygon 231 33 4.47665e-01 6.68e-07
polygon 232 19 1.34593e-01 2.01e-07
polygon 233 39 1.36327e+00 2.04e-06
polygon 234 137 1.55547e+02 2.32e-04
polygon 235 36 8.76479e+00 1.31e-05
polygon 236 79 3.08116e+01 4.60e-05
polygon 237 388 2.25271e+02 3.36e-04
polygon 238 316 7.78512e+01 1.16e-04
polygon 239 13 1.09564e-01 1.64e-07
polygon 240 18 3.49727e-01 5.22e-07
polygon 241 31 1.23017e+00 1.84e-06
polygon 242 16 6.55537e-01 9.79e-07
polygon 243 24 8.49487e-01 1.27e-06
polygon 244 30 2.54436e+00 3.80e-06
polygon 245 336 4.15806e+01 6.21e-05
polygon 246 330 1.69190e+02 2.53e-04
polygon 247 47 1.08035e+01 1.61e-05
polygon 248 39 4.94369e+00 7.38e-06
polygon 249 23 2.72438e+00 4.07e-06
polygon 250 33 5.70263e+00 8.52e-06
polygon 251 90 4.20329e+01 6.28e-05
polygon 252 28 1.35341e+00 2.02e-06
polygon 253 225 1.08816e+02 1.62e-04
polygon 254 33 9.16670e+00 1.37e-05
polygon 255 192 7.02655e+01 1.05e-04
polygon 256 49 1.49245e+01 2.23e-05
polygon 257 98 1.79076e+01 2.67e-05
polygon 258 6 6.37552e-01 9.52e-07
polygon 259 49 1.01233e+01 1.51e-05
polygon 260 141 3.43053e+01 5.12e-05
polygon 261 195 3.24345e+01 4.84e-05
polygon 262 51 3.38313e+00 5.05e-06
polygon 263 34 2.01400e+00 3.01e-06
polygon 264 13 2.50435e-01 3.74e-07
polygon 265 9 9.04824e-02 1.35e-07
polygon 266 34 4.61794e+00 6.90e-06
polygon 267 17 4.58200e-01 6.84e-07
polygon 268 15 2.74776e-01 4.10e-07
polygon 269 21 5.34978e-01 7.99e-07
polygon 270 19 4.55347e-01 6.80e-07
polygon 271 71 3.42557e+00 5.11e-06
polygon 272 24 1.32420e+00 1.98e-06
polygon 273 15 3.26247e-01 4.87e-07
polygon 274 39 8.65790e-01 1.29e-06
polygon 275 43 1.41627e+00 2.11e-06
polygon 276 24 7.52068e-01 1.12e-06
polygon 277 96 1.32101e+01 1.97e-05
polygon 278 38 1.18003e+00 1.76e-06
polygon 279 429 5.99087e+02 8.95e-04
polygon 280 13 1.74105e-01 2.60e-07
polygon 281 19 2.52336e-01 3.77e-07
polygon 282 16 3.11495e-01 4.65e-07
polygon 283 11 9.11047e-02 1.36e-07
polygon 284 12 2.13470e-01 3.19e-07
polygon 285 17 5.82663e-01 8.70e-07
polygon 286 56 2.60440e+01 3.89e-05
polygon 287 107 4.91389e+00 7.34e-06
polygon 288 51 2.79076e+00 4.17e-06
polygon 289 89 1.61156e+01 2.41e-05
polygon 290 28 1.30499e+00 1.95e-06
polygon 291 11 1.27616e-01 1.91e-07
polygon 292 34 2.54199e+00 3.80e-06
polygon 293 27 1.72476e+00 2.58e-06
polygon 294 37 2.01882e+00 3.01e-06
polygon 295 23 1.65571e+00 2.47e-06
polygon 296 33 3.05816e+00 4.57e-06
polygon 297 14 3.23153e-01 4.83e-07
polygon 298 91 1.51209e+01 2.26e-05
polygon 299 12 2.42901e-01 3.63e-07
polygon 300 11 1.37889e-01 2.06e-07
polygon 301 58 2.29751e+01 3.43e-05
polygon 302 48 5.10265e+00 7.62e-06
polygon 303 22 1.30706e+00 1.95e-06
polygon 304 15 3.49480e-01 5.22e-07
polygon 305 17 1.57570e+00 2.35e-06
polygon 306 34 3.68725e+00 5.51e-06
polygon 307 34 5.21904e+00 7.79e-06
polygon 308 24 5.42734e+00 8.10e-06
polygon 309 422 4.66497e+02 6.97e-04
polygon 310 142 2.98767e+01 4.46e-05
polygon 311 132 2.18707e+01 3.27e-05
polygon 312 19 5.88230e-01 8.78e-07
polygon 313 22 1.77611e+00 2.65e-06
polygon 314 (hole) 4 -6.42926e-10 -9.60e-16
polygon 315 40 4.09952e+00 6.12e-06
polygon 316 (hole) 3 -2.52076e-09 -3.76e-15
polygon 317 (hole) 24 -5.61620e-07 -8.39e-13
polygon 318 28 1.47685e+00 2.21e-06
polygon 319 (hole) 8 -8.93461e-08 -1.33e-13
polygon 320 (hole) 8 -1.34123e-07 -2.00e-13
polygon 321 67 9.99685e+00 1.49e-05
polygon 322 (hole) 8 -6.36345e-08 -9.50e-14
polygon 323 (hole) 6 -7.57978e-08 -1.13e-13
polygon 324 (hole) 3 -4.32022e-12 -6.45e-18
polygon 325 (hole) 4 -2.24224e-08 -3.35e-14
polygon 326 (hole) 5 -9.35663e-08 -1.40e-13
polygon 327 (hole) 13 -1.39251e-07 -2.08e-13
polygon 328 (hole) 7 -1.46533e-07 -2.19e-13
polygon 329 (hole) 6 -7.59071e-08 -1.13e-13
polygon 330 (hole) 7 -5.15915e-09 -7.70e-15
polygon 331 (hole) 10 -4.09465e-08 -6.11e-14
polygon 332 (hole) 21 -1.86624e-07 -2.79e-13
polygon 333 (hole) 3 -1.69684e-08 -2.53e-14
polygon 334 (hole) 8 -4.08363e-08 -6.10e-14
polygon 335 (hole) 3 -1.73631e-08 -2.59e-14
polygon 336 (hole) 5 -1.41820e-08 -2.12e-14
polygon 337 (hole) 10 -6.15226e-08 -9.19e-14
polygon 338 (hole) 8 -2.90782e-08 -4.34e-14
polygon 339 (hole) 8 -6.15579e-08 -9.19e-14
polygon 340 (hole) 4 -7.46371e-10 -1.11e-15
polygon 341 (hole) 4 -1.05432e-08 -1.57e-14
polygon 342 (hole) 4 -8.25148e-09 -1.23e-14
polygon 343 (hole) 3 -2.59991e-08 -3.88e-14
polygon 344 (hole) 3 -9.09758e-09 -1.36e-14
polygon 345 (hole) 3 -1.96869e-08 -2.94e-14
polygon 346 (hole) 3 -4.01481e-11 -5.99e-17
polygon 347 (hole) 5 -1.49421e-08 -2.23e-14
polygon 348 (hole) 4 -9.23861e-09 -1.38e-14
polygon 349 3 4.91147e-16 7.33e-22
polygon 350 (hole) 7 -5.73304e-08 -8.56e-14
polygon 351 (hole) 4 -4.26735e-08 -6.37e-14
polygon 352 (hole) 7 -1.69748e-08 -2.53e-14
polygon 353 (hole) 3 -1.05617e-08 -1.58e-14
polygon 354 25 4.82266e-01 7.20e-07
polygon 355 (hole) 6 -3.72405e-08 -5.56e-14
polygon 356 (hole) 7 -3.08332e-08 -4.60e-14
polygon 357 (hole) 4 -6.69857e-09 -1.00e-14
polygon 358 (hole) 4 -1.15937e-08 -1.73e-14
polygon 359 (hole) 4 -3.11046e-09 -4.64e-15
polygon 360 (hole) 12 -1.45689e-07 -2.18e-13
polygon 361 (hole) 4 -5.64800e-08 -8.43e-14
polygon 362 (hole) 3 -1.11188e-08 -1.66e-14
polygon 363 (hole) 4 -2.60006e-08 -3.88e-14
polygon 364 (hole) 28 -1.79082e-07 -2.67e-13
polygon 365 (hole) 4 -1.58213e-10 -2.36e-16
polygon 366 (hole) 4 -9.07739e-12 -1.36e-17
polygon 367 (hole) 4 -6.73612e-09 -1.01e-14
polygon 368 16 4.14093e-01 6.18e-07
polygon 369 (hole) 4 -5.36782e-09 -8.02e-15
polygon 370 (hole) 6 -1.51450e-08 -2.26e-14
polygon 371 (hole) 4 -4.38078e-09 -6.54e-15
polygon 372 (hole) 4 -1.36409e-08 -2.04e-14
polygon 373 (hole) 4 -3.84271e-08 -5.74e-14
polygon 374 (hole) 3 -7.56649e-09 -1.13e-14
polygon 375 (hole) 4 -8.04108e-09 -1.20e-14
polygon 376 (hole) 6 -2.53481e-08 -3.78e-14
polygon 377 (hole) 4 -1.15603e-08 -1.73e-14
polygon 378 (hole) 3 -2.16577e-09 -3.23e-15
polygon 379 (hole) 3 -2.87541e-09 -4.29e-15
polygon 380 (hole) 4 -5.42479e-08 -8.10e-14
polygon 381 (hole) 3 -4.75090e-09 -7.09e-15
polygon 382 (hole) 4 -3.45047e-08 -5.15e-14
polygon 383 (hole) 6 -7.20472e-09 -1.08e-14
polygon 384 (hole) 4 -2.01300e-08 -3.01e-14
polygon 385 (hole) 9 -7.99880e-08 -1.19e-13
polygon 386 (hole) 4 -1.30938e-09 -1.96e-15
polygon 387 (hole) 9 -1.92423e-08 -2.87e-14
polygon 388 (hole) 3 -2.44412e-13 -3.65e-19
polygon 389 (hole) 4 -2.73889e-08 -4.09e-14
polygon 390 (hole) 3 -1.15109e-08 -1.72e-14
polygon 391 (hole) 4 -8.77810e-08 -1.31e-13
polygon 392 26 3.03928e+00 4.54e-06
polygon 393 (hole) 6 -8.99139e-08 -1.34e-13
polygon 394 (hole) 7 -1.69870e-07 -2.54e-13
polygon 395 (hole) 11 -1.85737e-07 -2.77e-13
polygon 396 (hole) 32 -8.16597e-07 -1.22e-12
polygon 397 (hole) 7 -6.59191e-08 -9.84e-14
polygon 398 (hole) 3 -1.48932e-08 -2.22e-14
polygon 399 (hole) 13 -3.56626e-07 -5.33e-13
polygon 400 (hole) 3 -1.64295e-12 -2.45e-18
polygon 401 (hole) 3 -2.03106e-09 -3.03e-15
polygon 402 14 1.53563e-01 2.29e-07
polygon 403 (hole) 4 -1.79041e-08 -2.67e-14
polygon 404 (hole) 15 -2.31620e-07 -3.46e-13
polygon 405 (hole) 7 -3.46670e-08 -5.18e-14
polygon 406 (hole) 6 -1.90172e-08 -2.84e-14
polygon 407 (hole) 6 -3.15828e-08 -4.72e-14
polygon 408 (hole) 13 -1.40633e-08 -2.10e-14
polygon 409 (hole) 6 -2.76240e-08 -4.12e-14
polygon 410 (hole) 3 -3.74153e-08 -5.59e-14
polygon 411 (hole) 3 -1.16808e-08 -1.74e-14
polygon 412 (hole) 10 -4.10705e-09 -6.13e-15
polygon 413 (hole) 4 -6.52777e-09 -9.75e-15
polygon 414 (hole) 6 -5.03054e-08 -7.51e-14
polygon 415 (hole) 4 -2.54303e-09 -3.80e-15
polygon 416 (hole) 3 -1.76691e-10 -2.64e-16
polygon 417 (hole) 4 -7.90017e-08 -1.18e-13
polygon 418 (hole) 3 -1.07902e-09 -1.61e-15
polygon 419 (hole) 3 -8.74115e-10 -1.31e-15
polygon 420 (hole) 12 -6.41926e-08 -9.59e-14
polygon 421 (hole) 4 -4.36895e-09 -6.52e-15
polygon 422 (hole) 4 -1.85182e-08 -2.77e-14
polygon 423 (hole) 10 -1.21248e-07 -1.81e-13
polygon 424 (hole) 10 -5.59537e-08 -8.35e-14
polygon 425 (hole) 4 -2.20079e-08 -3.29e-14
polygon 426 19 5.08538e-01 7.59e-07
polygon 427 (hole) 4 -1.45573e-08 -2.17e-14
polygon 428 (hole) 4 -4.47181e-09 -6.68e-15
polygon 429 16 1.80565e-01 2.70e-07
polygon 430 (hole) 3 -9.03958e-14 -1.35e-19
polygon 431 26 9.75091e-01 1.46e-06
polygon 432 (hole) 4 -1.90718e-09 -2.85e-15
polygon 433 (hole) 5 -1.60264e-08 -2.39e-14
polygon 434 (hole) 4 -1.82919e-08 -2.73e-14
polygon 435 (hole) 3 -7.26182e-09 -1.08e-14
polygon 436 (hole) 4 -2.12571e-08 -3.17e-14
polygon 437 (hole) 4 -2.41922e-09 -3.61e-15
polygon 438 (hole) 3 -2.80771e-08 -4.19e-14
polygon 439 (hole) 7 -5.31797e-09 -7.94e-15
polygon 440 (hole) 16 -1.29582e-07 -1.93e-13
polygon 441 (hole) 4 -1.88745e-08 -2.82e-14
polygon 442 (hole) 3 -1.45605e-08 -2.17e-14
polygon 443 (hole) 3 -2.28810e-08 -3.42e-14
polygon 444 (hole) 4 -1.59810e-08 -2.39e-14
polygon 445 (hole) 5 -1.72617e-08 -2.58e-14
polygon 446 (hole) 3 -4.96859e-09 -7.42e-15
polygon 447 (hole) 14 -8.85135e-08 -1.32e-13
polygon 448 (hole) 4 -1.85807e-08 -2.77e-14
polygon 449 (hole) 4 -8.61250e-09 -1.29e-14
polygon 450 (hole) 4 -2.06032e-08 -3.08e-14
polygon 451 (hole) 3 -1.66645e-09 -2.49e-15
polygon 452 (hole) 7 -1.66808e-08 -2.49e-14
polygon 453 (hole) 4 -1.67544e-09 -2.50e-15
polygon 454 (hole) 12 -3.97402e-08 -5.93e-14
polygon 455 (hole) 4 -4.18178e-08 -6.24e-14
polygon 456 (hole) 3 -2.93135e-09 -4.38e-15
polygon 457 (hole) 4 -6.31983e-10 -9.44e-16
polygon 458 (hole) 3 -8.50840e-13 -1.27e-18
polygon 459 (hole) 7 -8.42406e-09 -1.26e-14
polygon 460 (hole) 3 -1.34855e-12 -2.01e-18
polygon 461 (hole) 6 -1.60663e-08 -2.40e-14
polygon 462 (hole) 4 -5.66656e-11 -8.46e-17
polygon 463 (hole) 5 -1.17874e-08 -1.76e-14
polygon 464 (hole) 3 -1.75180e-14 -2.62e-20
polygon 465 (hole) 5 -2.62273e-08 -3.92e-14
polygon 466 (hole) 7 -2.12912e-08 -3.18e-14
polygon 467 (hole) 11 -4.47347e-08 -6.68e-14
polygon 468 (hole) 6 -5.00912e-09 -7.48e-15
polygon 469 (hole) 5 -3.58831e-08 -5.36e-14
polygon 470 (hole) 3 -1.06446e-09 -1.59e-15
polygon 471 (hole) 4 -1.50038e-08 -2.24e-14
polygon 472 (hole) 3 -8.74496e-10 -1.31e-15
polygon 473 (hole) 7 -3.75991e-08 -5.61e-14
polygon 474 19 2.15295e-01 3.21e-07
polygon 475 (hole) 6 -3.71337e-08 -5.54e-14
polygon 476 (hole) 3 -4.33900e-09 -6.48e-15
polygon 477 (hole) 7 -8.65067e-08 -1.29e-13
polygon 478 (hole) 8 -7.76951e-08 -1.16e-13
polygon 479 (hole) 6 -9.01907e-08 -1.35e-13
polygon 480 (hole) 8 -1.10404e-07 -1.65e-13
polygon 481 (hole) 8 -7.38263e-08 -1.10e-13
polygon 482 (hole) 3 -1.71517e-08 -2.56e-14
polygon 483 16 3.50061e-01 5.23e-07
polygon 484 (hole) 3 -1.50971e-14 -2.25e-20
polygon 485 (hole) 3 -1.78834e-04 -2.67e-10
polygon 486 3 0.00000e+00 0.00e+00
polygon 487 (hole) 4 -4.21570e-08 -6.29e-14
polygon 488 (hole) 7 -1.74178e-07 -2.60e-13
polygon 489 14 2.11571e-01 3.16e-07
polygon 490 (hole) 3 -1.05578e-08 -1.58e-14
polygon 491 (hole) 4 -4.32122e-08 -6.45e-14
polygon 492 11 1.59536e-01 2.38e-07
polygon 493 (hole) 4 -3.72418e-08 -5.56e-14
polygon 494 (hole) 6 -2.87377e-08 -4.29e-14
polygon 495 (hole) 3 -3.87472e-09 -5.79e-15
polygon 496 (hole) 6 -1.05664e-08 -1.58e-14
polygon 497 (hole) 4 -7.04519e-09 -1.05e-14
polygon 498 111 7.35028e+01 1.10e-04
polygon 499 (hole) 3 -5.99345e-09 -8.95e-15
polygon 500 (hole) 3 -2.33491e-08 -3.49e-14
polygon 501 (hole) 5 -1.14503e-07 -1.71e-13
polygon 502 19 6.34649e-01 9.48e-07
polygon 503 (hole) 3 -2.26548e-09 -3.38e-15
polygon 504 (hole) 4 -5.35016e-08 -7.99e-14
polygon 505 10 3.74827e-02 5.60e-08
polygon 506 13 1.37100e-01 2.05e-07
polygon 507 31 4.47193e+00 6.68e-06
polygon 508 19 1.23742e+00 1.85e-06
polygon 509 20 1.79201e+00 2.68e-06
polygon 510 (hole) 3 -1.64235e-08 -2.45e-14
polygon 511 16 2.19464e-01 3.28e-07
polygon 512 11 7.94092e-02 1.19e-07
polygon 513 (hole) 3 -2.99572e-08 -4.47e-14
polygon 514 (hole) 3 -8.18639e-09 -1.22e-14
polygon 515 (hole) 3 -2.81986e-08 -4.21e-14
polygon 516 14 1.45779e-01 2.18e-07
polygon 517 31 6.24672e-01 9.33e-07
polygon 518 11 2.29288e-02 3.42e-08
polygon 519 18 2.10998e-01 3.15e-07
polygon 520 26 2.52571e+00 3.77e-06
polygon 521 16 5.77783e-01 8.63e-07
polygon 522 26 8.87985e-01 1.33e-06
polygon 523 13 1.10557e-01 1.65e-07
polygon 524 21 9.55338e-01 1.43e-06
polygon 525 32 2.18189e+00 3.26e-06
polygon 526 18 7.64830e-01 1.14e-06
polygon 527 34 1.71336e+00 2.56e-06
polygon 528 115 9.29594e-01 1.39e-06
polygon 529 68 2.06761e-01 3.09e-07
polygon 530 40 4.87851e+00 7.28e-06
polygon 531 10 7.43137e-02 1.11e-07
polygon 532 17 5.95020e-01 8.88e-07
polygon 533 21 4.02132e-01 6.00e-07
polygon 534 21 1.09635e+00 1.64e-06
polygon 535 14 1.15356e-01 1.72e-07
polygon 536 (hole) 18 -1.39516e-07 -2.08e-13
polygon 537 20 3.13610e-01 4.68e-07
polygon 538 275 3.92648e+02 5.86e-04
polygon 539 16 2.19612e-01 3.28e-07
polygon 540 13 2.64858e-01 3.95e-07
polygon 541 23 1.25549e+00 1.87e-06
polygon 542 9 1.79312e-01 2.68e-07
polygon 543 15 2.93132e-01 4.38e-07
polygon 544 15 5.29892e-01 7.91e-07
polygon 545 23 4.95351e-01 7.40e-07
polygon 546 (hole) 10 -2.23148e-07 -3.33e-13
polygon 547 14 4.41976e-01 6.60e-07
polygon 548 19 3.32318e-01 4.96e-07
polygon 549 (hole) 6 -6.70431e-08 -1.00e-13
polygon 550 16 4.15914e-01 6.21e-07
polygon 551 22 5.10526e-01 7.62e-07
polygon 552 48 1.87148e+01 2.79e-05
polygon 553 18 1.62882e+00 2.43e-06
polygon 554 (hole) 9 -1.96623e-07 -2.94e-13
polygon 555 71 2.43615e+01 3.64e-05
polygon 556 116 5.16444e+01 7.71e-05
polygon 557 19 1.54144e+00 2.30e-06
polygon 558 11 2.01013e-01 3.00e-07
polygon 559 (hole) 4 -1.87899e-08 -2.81e-14
polygon 560 20 2.31201e+00 3.45e-06
polygon 561 15 5.94674e-01 8.88e-07
polygon 562 218 1.56252e+02 2.33e-04
polygon 563 162 8.80247e+01 1.31e-04
polygon 564 16 3.69968e-01 5.52e-07
polygon 565 14 1.61451e-01 2.41e-07
polygon 566 18 5.54566e-01 8.28e-07
polygon 567 19 5.07078e-01 7.57e-07
polygon 568 (hole) 7 -1.85276e-07 -2.77e-13
polygon 569 57 6.86475e+00 1.03e-05
polygon 570 15 6.56018e-01 9.80e-07
polygon 571 21 6.65092e-01 9.93e-07
polygon 572 15 5.97784e-01 8.93e-07
polygon 573 37 2.56010e+00 3.82e-06
polygon 574 13 1.66996e-01 2.49e-07
polygon 575 20 8.18381e-01 1.22e-06
polygon 576 21 2.89704e+00 4.33e-06
polygon 577 15 4.08213e-01 6.10e-07
polygon 578 (hole) 3 -1.98618e-09 -2.97e-15
polygon 579 104 4.71547e+01 7.04e-05
polygon 580 28 3.80443e+00 5.68e-06
polygon 581 67 3.17343e+01 4.74e-05
polygon 582 730 9.20735e+02 1.37e-03
polygon 583 16 2.01576e-01 3.01e-07
polygon 584 12 1.41506e-01 2.11e-07
polygon 585 15 4.92125e-01 7.35e-07
polygon 586 18 1.58128e+00 2.36e-06
polygon 587 (hole) 7 -6.67005e-08 -9.96e-14
polygon 588 25 1.94322e+00 2.90e-06
polygon 589 12 2.51373e-01 3.75e-07
polygon 590 412 4.47936e+02 6.69e-04
polygon 591 22 7.57105e-01 1.13e-06
polygon 592 16 6.97434e-01 1.04e-06
polygon 593 (hole) 3 -2.86672e-08 -4.28e-14
polygon 594 39 4.45130e+00 6.65e-06
polygon 595 (hole) 4 -2.17002e-08 -3.24e-14
polygon 596 (hole) 3 -1.55741e-08 -2.33e-14
polygon 597 13 8.83786e-01 1.32e-06
polygon 598 10 2.06200e-01 3.08e-07
polygon 599 (hole) 3 -3.50430e-08 -5.23e-14
polygon 600 76 3.60497e+01 5.38e-05
polygon 601 (hole) 5 -5.41615e-08 -8.09e-14
polygon 602 13 4.56433e-01 6.82e-07
polygon 603 28 2.30613e+00 3.44e-06
polygon 604 (hole) 4 -1.59313e-08 -2.38e-14
polygon 605 27 2.50338e+00 3.74e-06
polygon 606 14 4.69238e-01 7.01e-07
polygon 607 10 1.69886e-01 2.54e-07
polygon 608 (hole) 8 -1.09145e-07 -1.63e-13
polygon 609 16 5.42822e-01 8.11e-07
polygon 610 10 1.17633e-01 1.76e-07
polygon 611 8 7.08579e-02 1.06e-07
polygon 612 144 7.85300e+01 1.17e-04
polygon 613 37 1.93477e+00 2.89e-06
polygon 614 (hole) 3 -1.47049e-09 -2.20e-15
polygon 615 25 5.59996e-01 8.36e-07
polygon 616 (hole) 8 -1.20535e-07 -1.80e-13
polygon 617 44 2.86031e+00 4.27e-06
polygon 618 26 1.26276e+00 1.89e-06
polygon 619 149 1.37840e+02 2.06e-04
polygon 620 18 4.84958e-01 7.24e-07
polygon 621 73 2.93195e+01 4.38e-05
polygon 622 (hole) 5 -1.99458e-08 -2.98e-14
polygon 623 45 6.87481e+00 1.03e-05
polygon 624 (hole) 4 -1.81549e-08 -2.71e-14
polygon 625 64 1.80880e+01 2.70e-05
polygon 626 12 8.76879e-01 1.31e-06
polygon 627 26 2.78381e+00 4.16e-06
polygon 628 35 8.72326e+00 1.30e-05
polygon 629 53 1.06237e+01 1.59e-05
polygon 630 26 5.40467e+00 8.07e-06
polygon 631 148 1.05037e+02 1.57e-04
polygon 632 21 8.32798e-01 1.24e-06
polygon 633 21 2.23023e+00 3.33e-06
polygon 634 8 6.30805e-01 9.42e-07
polygon 635 78 3.67603e+01 5.49e-05
polygon 636 18 1.64745e+00 2.46e-06
polygon 637 58 1.31747e+01 1.97e-05
polygon 638 94 1.16837e+01 1.74e-05
polygon 639 53 3.17801e+00 4.75e-06
polygon 640 137 8.46454e+00 1.26e-05
polygon 641 23 4.36507e-01 6.52e-07
polygon 642 27 9.72136e-01 1.45e-06
polygon 643 55 1.98854e+00 2.97e-06
polygon 644 48 1.02651e+01 1.53e-05
polygon 645 15 3.16153e-01 4.72e-07
polygon 646 22 8.93201e-01 1.33e-06
polygon 647 18 2.91538e-01 4.35e-07
polygon 648 59 3.70736e+00 5.54e-06
polygon 649 11 1.55205e-01 2.32e-07
polygon 650 14 1.87401e-01 2.80e-07
polygon 651 11 6.54128e-02 9.77e-08
polygon 652 11 8.64659e-02 1.29e-07
polygon 653 12 2.25129e-01 3.36e-07
polygon 654 30 2.72178e+00 4.06e-06
polygon 655 (hole) 3 -2.26252e-06 -3.38e-12
polygon 656 180 1.80314e+01 2.69e-05
polygon 657 14 4.01942e-01 6.00e-07
polygon 658 103 1.33467e+01 1.99e-05
polygon 659 14 2.07708e-01 3.10e-07
polygon 660 55 4.38623e+00 6.55e-06
polygon 661 41 4.69733e+00 7.01e-06
polygon 662 74 3.92633e+00 5.86e-06
polygon 663 31 2.79477e+00 4.17e-06
polygon 664 17 8.74607e-01 1.31e-06
polygon 665 26 4.50962e-01 6.73e-07
polygon 666 18 1.15065e+00 1.72e-06
polygon 667 164 1.48037e+01 2.21e-05
polygon 668 12 3.17290e-01 4.74e-07
polygon 669 69 2.59551e+01 3.88e-05
polygon 670 37 7.85963e-01 1.17e-06
polygon 671 14 3.37873e-01 5.05e-07
polygon 672 12 1.77980e-01 2.66e-07
polygon 673 21 6.46195e-01 9.65e-07
polygon 674 24 2.69998e+00 4.03e-06
polygon 675 30 1.00675e+00 1.50e-06
polygon 676 57 4.22747e+00 6.31e-06
polygon 677 23 8.63543e-01 1.29e-06
polygon 678 11 2.01925e-01 3.02e-07
polygon 679 95 2.88232e+01 4.30e-05
polygon 680 452 9.64208e+01 1.44e-04
polygon 681 77 1.72297e+01 2.57e-05
polygon 682 16 3.33522e-01 4.98e-07
polygon 683 19 6.82895e-01 1.02e-06
polygon 684 43 3.95062e+00 5.90e-06
polygon 685 29 1.00709e+00 1.50e-06
polygon 686 21 8.90014e-01 1.33e-06
polygon 687 22 7.93453e-01 1.18e-06
polygon 688 13 4.52813e-01 6.76e-07
polygon 689 30 1.49817e+00 2.24e-06
polygon 690 48 4.23996e+00 6.33e-06
polygon 691 37 1.20290e+00 1.80e-06
polygon 692 19 4.85805e-01 7.25e-07
polygon 693 46 2.18001e+00 3.26e-06
polygon 694 7 1.33723e-01 2.00e-07
polygon 695 54 1.82059e+01 2.72e-05
polygon 696 10 8.57866e-01 1.28e-06
polygon 697 14 1.98445e-01 2.96e-07
polygon 698 19 4.68357e-01 6.99e-07
polygon 699 52 8.57765e+00 1.28e-05
polygon 700 73 7.91998e+00 1.18e-05
polygon 701 23 8.19561e+00 1.22e-05
polygon 702 10 1.06594e-01 1.59e-07
polygon 703 169 1.39462e+01 2.08e-05
polygon 704 17 2.55915e-01 3.82e-07
polygon 705 12 1.81516e-01 2.71e-07
polygon 706 16 3.68509e-01 5.50e-07
polygon 707 12 1.24954e-01 1.87e-07
polygon 708 1012 4.61886e+02 6.90e-04
polygon 709 15 5.42647e-01 8.10e-07
polygon 710 12 5.53088e-01 8.26e-07
polygon 711 24 1.06455e+00 1.59e-06
polygon 712 14 3.17038e-01 4.73e-07
polygon 713 17 1.25559e+00 1.87e-06
polygon 714 16 1.92407e+00 2.87e-06
polygon 715 80 1.86306e+01 2.78e-05
polygon 716 11 6.83880e-02 1.02e-07
polygon 717 166 4.21470e+01 6.29e-05
polygon 718 22 1.36174e+00 2.03e-06
polygon 719 31 6.07363e-01 9.07e-07
polygon 720 10 1.27067e-01 1.90e-07
polygon 721 164 5.02951e+01 7.51e-05
polygon 722 17 1.20798e+00 1.80e-06
polygon 723 59 2.41720e+00 3.61e-06
polygon 724 12 2.52952e-01 3.78e-07
polygon 725 14 2.44407e-01 3.65e-07
polygon 726 43 4.89177e+00 7.30e-06
polygon 727 40 3.86642e+00 5.77e-06
polygon 728 23 4.59203e-01 6.86e-07
polygon 729 14 1.79442e-01 2.68e-07
polygon 730 13 1.61191e-01 2.41e-07
polygon 731 43 1.33365e+00 1.99e-06
polygon 732 75 1.20247e+01 1.80e-05
polygon 733 17 9.96466e-02 1.49e-07
polygon 734 22 1.60660e+00 2.40e-06
polygon 735 730 1.16983e+02 1.75e-04
polygon 736 129 7.87539e+00 1.18e-05
polygon 737 22 4.66256e-01 6.96e-07
polygon 738 19 5.44519e-01 8.13e-07
polygon 739 16 3.73911e-01 5.58e-07
polygon 740 88 3.47678e+01 5.19e-05
polygon 741 46 3.23711e+00 4.83e-06
polygon 742 47 9.17482e-01 1.37e-06
polygon 743 14 4.44869e-01 6.64e-07
polygon 744 43 1.29527e+00 1.93e-06
polygon 745 65 3.20974e+00 4.79e-06
polygon 746 306 1.83695e+01 2.74e-05
polygon 747 17 4.25486e-01 6.35e-07
polygon 748 29 1.08467e+00 1.62e-06
polygon 749 26 1.16602e+00 1.74e-06
polygon 750 135 4.90655e+00 7.33e-06
polygon 751 23 1.87723e+00 2.80e-06
polygon 752 21 7.02338e-01 1.05e-06
polygon 753 30 4.99970e+00 7.47e-06
polygon 754 21 7.71571e-01 1.15e-06
polygon 755 13 1.47832e-01 2.21e-07
polygon 756 26 3.17307e+00 4.74e-06
polygon 757 17 1.03642e+00 1.55e-06
polygon 758 100 8.74537e+00 1.31e-05
polygon 759 11 1.13950e-01 1.70e-07
polygon 760 26 8.27840e-01 1.24e-06
polygon 761 17 4.10882e-01 6.14e-07
polygon 762 13 1.49250e-01 2.23e-07
polygon 763 35 1.76408e+00 2.63e-06
polygon 764 15 2.03390e-01 3.04e-07
polygon 765 18 1.70982e-01 2.55e-07
polygon 766 32 7.98599e+00 1.19e-05
polygon 767 27 3.20839e-01 4.79e-07
polygon 768 539 1.17139e+02 1.75e-04
polygon 769 51 1.57614e+00 2.35e-06
polygon 770 19 3.60772e-01 5.39e-07
polygon 771 58 2.10060e+00 3.14e-06
polygon 772 13 1.70224e-01 2.54e-07
polygon 773 11 6.11736e-01 9.13e-07
polygon 774 35 5.26152e+00 7.86e-06
polygon 775 16 1.99371e-01 2.98e-07
polygon 776 14 4.48979e-01 6.70e-07
polygon 777 25 1.92913e+00 2.88e-06
polygon 778 16 4.26157e-01 6.36e-07
polygon 779 46 1.88229e+00 2.81e-06
polygon 780 17 2.29456e-01 3.43e-07
polygon 781 36 4.01418e+00 5.99e-06
polygon 782 57 1.77153e+01 2.65e-05
polygon 783 20 5.72298e-01 8.55e-07
polygon 784 14 2.46782e-01 3.68e-07
polygon 785 29 1.01300e+00 1.51e-06
polygon 786 90 8.54955e+00 1.28e-05
polygon 787 21 3.99098e-01 5.96e-07
polygon 788 12 2.81650e-01 4.21e-07
polygon 789 52 3.19337e+00 4.77e-06
polygon 790 13 3.40300e-01 5.08e-07
polygon 791 29 9.08533e-01 1.36e-06
polygon 792 20 3.07393e-01 4.59e-07
polygon 793 19 1.15906e+00 1.73e-06
polygon 794 17 1.16182e+00 1.73e-06
polygon 795 21 2.60877e+00 3.90e-06
polygon 796 15 1.09001e-01 1.63e-07
polygon 797 37 1.31921e+00 1.97e-06
polygon 798 53 1.46854e+00 2.19e-06
polygon 799 59 8.92822e+00 1.33e-05
polygon 800 9 7.44981e-02 1.11e-07
polygon 801 14 2.41462e-01 3.61e-07
polygon 802 96 7.25940e+00 1.08e-05
polygon 803 11 1.06055e-01 1.58e-07
polygon 804 49 1.87834e+00 2.80e-06
polygon 805 23 6.08310e-01 9.08e-07
polygon 806 50 6.76488e-01 1.01e-06
polygon 807 22 4.83089e-01 7.21e-07
polygon 808 17 1.17278e-01 1.75e-07
polygon 809 13 8.98786e-01 1.34e-06
polygon 810 43 1.01757e+00 1.52e-06
polygon 811 52 1.68377e+00 2.51e-06
polygon 812 348 2.50314e+02 3.74e-04
polygon 813 43 1.29120e+00 1.93e-06
polygon 814 71 2.02836e+00 3.03e-06
polygon 815 (hole) 4 -4.12791e-06 -6.16e-12
polygon 816 141 4.53240e+00 6.77e-06
polygon 817 52 3.62008e+00 5.41e-06
polygon 818 20 7.69539e-01 1.15e-06
polygon 819 88 8.88904e+00 1.33e-05
polygon 820 12 1.44668e-01 2.16e-07
polygon 821 46 2.28026e+00 3.40e-06
polygon 822 39 4.79165e+00 7.15e-06
polygon 823 53 2.86736e+00 4.28e-06
polygon 824 38 3.22508e+00 4.82e-06
polygon 825 18 4.44863e-01 6.64e-07
polygon 826 45 1.45134e+00 2.17e-06
polygon 827 86 2.56400e+01 3.83e-05
polygon 828 158 1.74510e+01 2.61e-05
polygon 829 14 2.50383e-01 3.74e-07
polygon 830 120 7.17019e+00 1.07e-05
polygon 831 146 2.28448e+01 3.41e-05
polygon 832 143 2.34472e+01 3.50e-05
polygon 833 20 3.47415e-01 5.19e-07
polygon 834 (hole) 4 -4.70080e-12 -7.02e-18
polygon 835 84 3.86387e+00 5.77e-06
polygon 836 68 1.10850e+01 1.66e-05
polygon 837 863 7.47702e+01 1.12e-04
polygon 838 53 7.54742e+00 1.13e-05
polygon 839 (hole) 5 -5.07266e-12 -7.57e-18
polygon 840 64 2.11016e+00 3.15e-06
polygon 841 103 1.47525e+01 2.20e-05
polygon 842 26 3.40345e+00 5.08e-06
polygon 843 (hole) 6 -3.63322e-12 -5.43e-18
polygon 844 (hole) 8 -7.62695e-12 -1.14e-17
polygon 845 103 8.40445e+00 1.25e-05
polygon 846 23 4.92041e-01 7.35e-07
polygon 847 27 5.53625e-01 8.27e-07
polygon 848 60 3.56102e+00 5.32e-06
polygon 849 151 1.51314e+01 2.26e-05
polygon 850 (hole) 3 -7.39478e-13 -1.10e-18
polygon 851 708 4.18129e+02 6.24e-04
polygon 852 120 1.32857e+01 1.98e-05
polygon 853 856 2.53898e+02 3.79e-04
polygon 854 (hole) 3 -1.51167e-12 -2.26e-18
polygon 855 (hole) 3 -4.16327e-13 -6.22e-19
polygon 856 (hole) 3 -1.81968e-13 -2.72e-19
polygon 857 (hole) 5 -8.69542e-12 -1.30e-17
polygon 858 13 2.47298e-01 3.69e-07
polygon 859 (hole) 5 -7.23071e-12 -1.08e-17
polygon 860 16 7.35754e-01 1.10e-06
polygon 861 (hole) 4 -5.36614e-12 -8.01e-18
polygon 862 21 5.31633e-01 7.94e-07
polygon 863 (hole) 11 -4.03802e-12 -6.03e-18
polygon 864 (hole) 8 -8.21058e-12 -1.23e-17
polygon 865 (hole) 3 -8.94235e-13 -1.34e-18
polygon 866 (hole) 5 -8.77891e-13 -1.31e-18
polygon 867 72 1.03185e+01 1.54e-05
polygon 868 (hole) 4 -1.46473e-12 -2.19e-18
polygon 869 (hole) 3 -9.57738e-14 -1.43e-19
polygon 870 38 1.34734e+00 2.01e-06
polygon 871 8 4.09489e-01 6.11e-07
polygon 872 (hole) 4 -2.72557e-12 -4.07e-18
polygon 873 17 2.46995e-01 3.69e-07
polygon 874 (hole) 3 -4.56587e-13 -6.82e-19
polygon 875 14 1.42285e-01 2.12e-07
polygon 876 (hole) 4 -1.74558e-12 -2.61e-18
polygon 877 (hole) 6 -5.91253e-12 -8.83e-18
polygon 878 (hole) 4 -1.34510e-13 -2.01e-19
polygon 879 (hole) 3 -6.12765e-17 -9.15e-23
polygon 880 9 1.56346e-01 2.33e-07
polygon 881 4 2.90958e-17 4.34e-23
polygon 882 42 5.49581e-01 8.21e-07
polygon 883 (hole) 5 -2.12744e-12 -3.18e-18
polygon 884 (hole) 3 -5.23582e-12 -7.82e-18
polygon 885 (hole) 3 -2.43443e-13 -3.64e-19
polygon 886 (hole) 3 -8.02710e-13 -1.20e-18
polygon 887 (hole) 6 -1.09890e-11 -1.64e-17
polygon 888 (hole) 4 -1.53858e-13 -2.30e-19
polygon 889 59 1.15813e+00 1.73e-06
polygon 890 (hole) 4 -2.30355e-12 -3.44e-18
polygon 891 12 1.16845e-01 1.74e-07
polygon 892 11 6.00122e-02 8.96e-08
polygon 893 14 1.63465e-01 2.44e-07
polygon 894 8 4.73674e-02 7.07e-08
polygon 895 23 6.69330e-01 9.99e-07
polygon 896 (hole) 3 -1.13772e-13 -1.70e-19
polygon 897 144 6.48346e+00 9.68e-06
polygon 898 74 4.79923e+00 7.17e-06
polygon 899 116 1.02706e+01 1.53e-05
polygon 900 16 2.74005e-01 4.09e-07
polygon 901 97 4.72599e+00 7.06e-06
polygon 902 18 6.15583e-01 9.19e-07
polygon 903 146 1.38979e+01 2.08e-05
polygon 904 (hole) 5 -4.08205e-12 -6.10e-18
polygon 905 32 1.94516e+00 2.90e-06
polygon 906 26 9.81356e-01 1.47e-06
polygon 907 (hole) 3 -2.01834e-12 -3.01e-18
polygon 908 (hole) 4 -9.94985e-13 -1.49e-18
polygon 909 82 4.39797e+00 6.57e-06
polygon 910 33 9.36113e-01 1.40e-06
polygon 911 (hole) 3 -1.40697e-13 -2.10e-19
polygon 912 47 1.51844e+00 2.27e-06
polygon 913 (hole) 3 -6.88122e-12 -1.03e-17
polygon 914 196 2.32689e+01 3.47e-05
polygon 915 (hole) 3 -2.40592e-12 -3.59e-18
polygon 916 (hole) 5 -6.19998e-12 -9.26e-18
polygon 917 53 1.30069e+00 1.94e-06
polygon 918 70 1.00802e+01 1.51e-05
polygon 919 (hole) 3 -3.17405e-13 -4.74e-19
polygon 920 10 5.33870e-01 7.97e-07
polygon 921 211 2.73026e+01 4.08e-05
polygon 922 41 1.44141e+00 2.15e-06
polygon 923 9 5.20220e-01 7.77e-07
polygon 924 316 6.66808e+01 9.96e-05
polygon 925 17 7.48505e-01 1.12e-06
polygon 926 15 1.49799e-01 2.24e-07
polygon 927 24 4.99500e-01 7.46e-07
polygon 928 12 2.24463e-01 3.35e-07
polygon 929 100 1.06028e+01 1.58e-05
polygon 930 11 1.24267e-01 1.86e-07
polygon 931 143 1.89581e+01 2.83e-05
polygon 932 60 7.07802e+00 1.06e-05
polygon 933 124 1.23476e+01 1.84e-05
polygon 934 34 5.67625e+00 8.48e-06
polygon 935 797 1.86955e+02 2.79e-04
polygon 936 15 7.15891e-01 1.07e-06
polygon 937 20 1.61624e+00 2.41e-06
polygon 938 13 1.73955e-01 2.60e-07
polygon 939 60 2.42590e+00 3.62e-06
polygon 940 20 3.86365e-01 5.77e-07
polygon 941 213 3.69421e+01 5.52e-05
polygon 942 15 4.55574e-01 6.80e-07
polygon 943 35 1.56058e+00 2.33e-06
polygon 944 32 1.42419e+00 2.13e-06
polygon 945 68 2.30711e+00 3.44e-06
polygon 946 240 5.93935e+01 8.87e-05
polygon 947 145 1.58875e+01 2.37e-05
polygon 948 39 3.67239e+00 5.48e-06
polygon 949 54 3.84847e+00 5.75e-06
polygon 950 32 1.14073e+00 1.70e-06
polygon 951 18 3.58549e-01 5.35e-07
polygon 952 25 8.06470e-01 1.20e-06
polygon 953 15 2.67818e-01 4.00e-07
polygon 954 26 1.43681e+00 2.15e-06
polygon 955 18 3.10061e-01 4.63e-07
polygon 956 18 4.71644e-01 7.04e-07
polygon 957 16 5.04146e-01 7.53e-07
polygon 958 72 7.33720e+00 1.10e-05
polygon 959 27 1.37772e+00 2.06e-06
polygon 960 15 3.17217e-01 4.74e-07
polygon 961 37 1.17498e+00 1.75e-06
polygon 962 21 7.05388e-01 1.05e-06
polygon 963 216 2.05399e+01 3.07e-05
polygon 964 16 2.40093e-01 3.59e-07
polygon 965 29 1.71282e+00 2.56e-06
polygon 966 30 1.28695e+00 1.92e-06
polygon 967 36 2.29670e+00 3.43e-06
polygon 968 129 1.97698e+01 2.95e-05
polygon 969 65 2.65969e+00 3.97e-06
polygon 970 283 3.30575e+01 4.94e-05
polygon 971 28 7.06271e-01 1.05e-06
polygon 972 28 5.50158e-01 8.21e-07
polygon 973 49 2.25015e+00 3.36e-06
polygon 974 26 1.24280e+00 1.86e-06
polygon 975 28 9.02794e-01 1.35e-06
polygon 976 18 4.84064e-01 7.23e-07
polygon 977 47 2.33442e+00 3.49e-06
polygon 978 17 2.23835e-01 3.34e-07
polygon 979 34 7.32504e-01 1.09e-06
polygon 980 41 7.15538e-01 1.07e-06
polygon 981 21 5.81687e-01 8.69e-07
polygon 982 63 3.83819e+00 5.73e-06
polygon 983 27 1.33192e+00 1.99e-06
polygon 984 46 2.95028e+00 4.41e-06
polygon 985 10 1.34210e-01 2.00e-07
polygon 986 16 2.38274e-01 3.56e-07
polygon 987 18 4.18536e-01 6.25e-07
polygon 988 46 1.26584e+00 1.89e-06
polygon 989 14 2.14679e-01 3.21e-07
polygon 990 76 3.63371e+00 5.43e-06
polygon 991 339 4.44685e+01 6.64e-05
polygon 992 109 6.29786e+00 9.40e-06
polygon 993 14 2.89570e-01 4.32e-07
polygon 994 37 9.07704e-01 1.36e-06
polygon 995 68 3.86104e+00 5.77e-06
polygon 996 183 1.45458e+01 2.17e-05
polygon 997 44 1.64355e+00 2.45e-06
polygon 998 26 1.21667e+00 1.82e-06
polygon 999 13 1.95710e-01 2.92e-07
polygon 1000 38 1.44402e+00 2.16e-06
polygon 1001 14 3.29394e-01 4.92e-07
polygon 1002 14 1.79828e-01 2.69e-07
polygon 1003 44 2.11118e+00 3.15e-06
polygon 1004 643 1.79215e+02 2.68e-04
polygon 1005 24 1.10979e+00 1.66e-06
polygon 1006 30 1.57106e+00 2.35e-06
polygon 1007 46 3.61716e+00 5.40e-06
polygon 1008 93 7.56340e-01 1.13e-06
polygon 1009 57 3.07732e+00 4.59e-06
polygon 1010 18 2.16224e-01 3.23e-07
polygon 1011 17 5.12388e-01 7.65e-07
polygon 1012 9 7.44780e-02 1.11e-07
polygon 1013 57 8.09122e+00 1.21e-05
polygon 1014 16 4.01323e-01 5.99e-07
polygon 1015 118 2.90304e+01 4.33e-05
polygon 1016 22 4.44693e-01 6.64e-07
polygon 1017 43 2.28763e+00 3.42e-06
polygon 1018 27 1.15444e+00 1.72e-06
polygon 1019 25 1.21695e+00 1.82e-06
polygon 1020 53 4.49228e+00 6.71e-06
polygon 1021 (hole) 8 -5.14959e-08 -7.69e-14
polygon 1022 (hole) 3 -1.87101e-12 -2.79e-18
polygon 1023 (hole) 31 -2.12239e-07 -3.17e-13
polygon 1024 (hole) 6 -4.99400e-08 -7.46e-14
polygon 1025 (hole) 8 -1.14003e-07 -1.70e-13
polygon 1026 (hole) 4 -5.05161e-08 -7.54e-14
polygon 1027 (hole) 4 -6.87492e-09 -1.03e-14
polygon 1028 (hole) 6 -5.81106e-12 -8.68e-18
polygon 1029 (hole) 7 -1.10596e-07 -1.65e-13
polygon 1030 (hole) 14 -3.52603e-07 -5.26e-13
polygon 1031 (hole) 14 -1.57740e-07 -2.36e-13
polygon 1032 (hole) 6 -6.42551e-08 -9.59e-14
polygon 1033 (hole) 4 -3.16857e-08 -4.73e-14
polygon 1034 (hole) 8 -7.25950e-08 -1.08e-13
polygon 1035 (hole) 4 -1.10823e-08 -1.65e-14
polygon 1036 (hole) 8 -1.54074e-07 -2.30e-13
polygon 1037 (hole) 4 -5.15531e-09 -7.70e-15
polygon 1038 (hole) 7 -4.08919e-08 -6.11e-14
polygon 1039 (hole) 9 -6.45354e-08 -9.64e-14
polygon 1040 (hole) 13 -6.64876e-08 -9.93e-14
polygon 1041 (hole) 5 -5.31928e-08 -7.94e-14
polygon 1042 (hole) 6 -4.67379e-08 -6.98e-14
polygon 1043 (hole) 4 -2.58840e-08 -3.86e-14
polygon 1044 (hole) 10 -9.35297e-08 -1.40e-13
polygon 1045 (hole) 4 -1.76462e-08 -2.63e-14
polygon 1046 (hole) 4 -5.14884e-08 -7.69e-14
polygon 1047 (hole) 3 -1.43941e-09 -2.15e-15
polygon 1048 (hole) 8 -2.11864e-08 -3.16e-14
polygon 1049 (hole) 4 -3.97950e-09 -5.94e-15
polygon 1050 (hole) 3 -2.72931e-08 -4.08e-14
polygon 1051 (hole) 4 -9.91931e-08 -1.48e-13
polygon 1052 (hole) 9 -4.30090e-08 -6.42e-14
polygon 1053 (hole) 3 -2.69442e-08 -4.02e-14
polygon 1054 (hole) 6 -3.65336e-08 -5.46e-14
polygon 1055 (hole) 10 -5.24779e-08 -7.84e-14
polygon 1056 (hole) 12 -5.88191e-08 -8.78e-14
polygon 1057 (hole) 8 -3.34285e-08 -4.99e-14
polygon 1058 (hole) 18 -1.75177e-07 -2.62e-13
polygon 1059 (hole) 6 -6.63687e-08 -9.91e-14
polygon 1060 (hole) 3 -1.24000e-08 -1.85e-14
polygon 1061 (hole) 10 -1.14560e-07 -1.71e-13
polygon 1062 (hole) 5 -5.13944e-08 -7.67e-14
polygon 1063 (hole) 10 -8.30360e-08 -1.24e-13
polygon 1064 (hole) 4 -1.89513e-08 -2.83e-14
polygon 1065 (hole) 6 -1.07333e-07 -1.60e-13
polygon 1066 (hole) 3 -5.80811e-09 -8.67e-15
polygon 1067 79 1.47390e+01 2.20e-05
polygon 1068 54 9.35305e+00 1.40e-05
polygon 1069 (hole) 3 -1.94488e-09 -2.90e-15
polygon 1070 (hole) 10 -6.90517e-08 -1.03e-13
polygon 1071 (hole) 4 -1.17368e-08 -1.75e-14
polygon 1072 (hole) 9 -5.79029e-08 -8.65e-14
polygon 1073 (hole) 3 -4.78872e-08 -7.15e-14
polygon 1074 (hole) 3 -5.90058e-09 -8.81e-15
polygon 1075 (hole) 5 -6.11108e-08 -9.12e-14
polygon 1076 37606 6.60254e+05 9.86e-01
polygon 1077 (hole) 5 -2.16281e-12 -3.23e-18
enclosing rectangle: [-210.0086, 724.6476] x [1072.0263, 3158.4671] km
(934.7 x 2086 km)
Window area = 669714 square km
Unit of length: 1 km
Fraction of frame area: 0.343The ppp object outputted from combining both the point and polygon feature results in the boundary of Myanmar outlining the plot of conflict events as shown.
# Plot each masked ppp object
par(mfrow = c(2,3), mar = c(0,0,1,0))
for (quarter in names(masked_ppp_list_km)) {
plot(masked_ppp_list_km[[quarter]], main = paste("Year Quarter:", quarter))
}
Let’s also create combined ppp and owin objects by the districts I am interested in this study. We’ll prepare the objects for the districts: Yinmarbin, Shwebo, Pakokku and Mandalay.
# Prepare Dataset for Yinmarbin District
conflict_yinmarbin = filter(conflict_data_sf, DT == "Yinmarbin")
boundary_yinmarbin <- filter(boundary_sf, DT == "Yinmarbin")
# Create a combined ppp and owin object
yinmarbin_owin <- as.owin(boundary_yinmarbin)
ppp_obj <- as.ppp(conflict_yinmarbin$geometry)
ppp_obj <- as.ppp(st_geometry(conflict_yinmarbin))
# Handle duplicates
ppp_obj <- rjitter(ppp_obj, retry=TRUE, nsim=1, drop=TRUE)
yinmarbin_ppp_owin <- ppp_obj[yinmarbin_owin]
yinmarbin_ppp_owin <- rescale(yinmarbin_ppp_owin, 1000, "km")# Prepare Dataset for Shwebo District
conflict_shwebo = filter(conflict_data_sf, DT == "Shwebo")
boundary_shwebo <- filter(boundary_sf, DT == "Shwebo")
# Create a combined ppp and owin object
shwebo_owin <- as.owin(boundary_shwebo)
ppp_obj <- as.ppp(conflict_shwebo$geometry)
# Handle duplicates
ppp_obj <- rjitter(ppp_obj, retry=TRUE, nsim=1, drop=TRUE)
shwebo_ppp_owin <- ppp_obj[shwebo_owin]
shwebo_ppp_owin <- rescale(shwebo_ppp_owin, 1000, "km")# Prepare Dataset for Pakokku District
conflict_pakokku = filter(conflict_data_sf, DT == "Pakokku")
boundary_pakokku <- filter(boundary_sf, DT == "Pakokku")
# Create a combined ppp and owin object
pakokku_owin <- as.owin(boundary_pakokku)
ppp_obj <- as.ppp(conflict_pakokku$geometry)
# Handle duplicates
ppp_obj <- rjitter(ppp_obj, retry=TRUE, nsim=1, drop=TRUE)
pakokku_ppp_owin <- ppp_obj[pakokku_owin]
pakokku_ppp_owin <- rescale(pakokku_ppp_owin, 1000, "km")# Prepare Dataset for Mandalay District
conflict_mandalay = filter(conflict_data_sf, DT == "Mandalay")
boundary_mandalay <- filter(boundary_sf, DT == "Mandalay")
# Create a combined ppp and owin object
mandalay_owin <- as.owin(boundary_mandalay)
ppp_obj <- as.ppp(conflict_mandalay$geometry)
# Handle duplicates
ppp_obj <- rjitter(ppp_obj, retry=TRUE, nsim=1, drop=TRUE)
mandalay_ppp_owin <- ppp_obj[mandalay_owin]
mandalay_ppp_owin <- rescale(mandalay_ppp_owin, 1000, "km")par(mfrow = c(2,2), mar = c(0,0,1,0))
plot(yinmarbin_ppp_owin, main = paste("Yinmarbin District"))
plot(shwebo_ppp_owin, main = paste("Shwebo District"))
plot(pakokku_ppp_owin, main = paste("Pakokku District"))
plot(mandalay_ppp_owin, main = paste("Mandalay District"))
Observations
Here I plot the districts with the highest concentration of armed conflicts. We can observe localised clustering of armed conflicts in certain towns of each district (darker spots on the map) which also indicates high frequency of conflicts occurring in these areas. On the other hand, there are also parts of each district with smaller frequency of conflicts but clusters of conflicts formed.
4. Exploratory Data Analysis
4.1 Identifying Districts with Highest Proportion of Conflicts
It’ll also be interesting to find out specific districts with the highest concentration of armed conflicts. I will first calculate the total occurrences of conflict events per district and add the column to boundary_sf.
Count number of conflicts by districts
conflict_count <- conflict_data_sf %>%
group_by(DT) %>%
summarise(total_count_DT = n()) %>%
st_drop_geometry() %>%
select(DT, total_count_DT)
# Perform the join
boundary_sf <- boundary_sf %>%
left_join(conflict_count, by = "DT")Next, let’s calculate the proportion of total conflicts and add it as a column into the boundary_sf dataset as proportion_DT.
# Create new 'proportion_DT' column
boundary_sf <- boundary_sf %>%
mutate(proportion_DT = total_count_DT / sum(total_count_DT))
head(boundary_sf[c('DT','total_count_DT','proportion_DT')])Simple feature collection with 6 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -14915.04 ymin: 1736124 xmax: 187961.7 ymax: 2051144
Projected CRS: WGS 84 / UTM zone 47N
DT total_count_DT proportion_DT geometry
1 Hinthada 160 NA MULTIPOLYGON (((90859.89 20...
2 Labutta 51 NA MULTIPOLYGON (((75991.51 17...
3 Maubin 118 NA MULTIPOLYGON (((115559 1928...
4 Myaungmya 59 NA MULTIPOLYGON (((39919.39 18...
5 Pathein 334 NA MULTIPOLYGON (((-6302.348 1...
6 Pyapon 131 NA MULTIPOLYGON (((93411.72 17...At a quick glance, we can see that central and southern parts of Myanmar have the highest proportions of armed conflict events occurring.
Set up the points map
districts_choropleth <-
tm_shape(boundary_sf) +
tm_fill("proportion_DT",
n=10,
title="Proportion",
style="equal",
palette="Blues") +
tm_borders(lwd=0.2,
alpha=1) +
tm_text(text = "DT",
size = 0.2,
col = "black",
fontface = "bold") +
tm_layout(main.title = "Distribution of Conflict Points Across Districts",
legend.outside=FALSE,
main.title.size=1)# Plot the map
tmap_mode("plot")
tmap_arrange(districts_choropleth)
More specifically, we can observe that conflict hotspots are mainly found in the districts of Yinmarbin, Shwebo, Pakokku and Mandalay which lies in the central regions of Myanmar.
4.2 Identifying States with Highest Proportion of Conflicts
Instead, let us also explore the top 10 states with the highest proportions of armed conflict events.
Count number of conflicts by states
conflict_count <- conflict_data_sf %>%
group_by(ST) %>%
summarise(total_count_ST = n()) %>%
st_drop_geometry() %>%
select(ST, total_count_ST)
# Perform the join
boundary_sf <- boundary_sf %>%
left_join(conflict_count, by = "ST")Likewise, I’ll add a new column called proportion_ST to represent the proportion based on each Myanmar state.
boundary_sf <- boundary_sf %>%
mutate(proportion_ST = total_count_ST / sum(total_count_ST))
head(boundary_sf[c('ST','total_count_ST','proportion_ST')])Simple feature collection with 6 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -14915.04 ymin: 1736124 xmax: 187961.7 ymax: 2051144
Projected CRS: WGS 84 / UTM zone 47N
ST total_count_ST proportion_ST geometry
1 Ayeyarwady 853 NA MULTIPOLYGON (((90859.89 20...
2 Ayeyarwady 853 NA MULTIPOLYGON (((75991.51 17...
3 Ayeyarwady 853 NA MULTIPOLYGON (((115559 1928...
4 Ayeyarwady 853 NA MULTIPOLYGON (((39919.39 18...
5 Ayeyarwady 853 NA MULTIPOLYGON (((-6302.348 1...
6 Ayeyarwady 853 NA MULTIPOLYGON (((93411.72 17...At a quick glance, we can see that central and southern parts of Myanmar have the highest proportions of armed conflict events occurring, namely in Sagaing, Mandalay, Magway and Yangon states as indicated in the map below (darkest shade of blue).
Create the points map
states_choropleth <-
tm_shape(boundary_sf) +
tm_fill("proportion_ST",
n=10,
title="Proportion",
style="equal",
palette="Blues") +
tm_borders(lwd=0.2,
alpha=1) +
tm_text(text = "ST",
size = 0.2,
col = "black",
fontface = "bold") +
tm_layout(main.title = "Distribution of Conflict Points Across States",
legend.outside=FALSE,
main.title.size=1)# Plot the map
tmap_mode("plot")
tmap_arrange(states_choropleth)
4.3 Standard Distances of Top 4 Districts
Next, I’ve plot the standard distances for the four districts that I’m most interested in. The codes here will calculate the average location of the conflict points within the district (i.e. the mean centre) and measure the dispersion of points around the mean centre.
The circle drawn around the mean centre shows us the standard distance to which we can observe that Mandalay has the smallest standard distance and hence, has the most conflict points closely clustered around the mean centre.
Find hot spot of Yinmarbin District
# Define the districts
districts <- c("Yinmarbin", "Shwebo", "Pakokku", "Mandalay")
for (district in districts) {
# Filter the conflict data for the current district
conflict_district <- filter(conflict_data_sf, DT == district)
boundary_district <- filter(boundary_sf, DT == district)
# Create a combined ppp and owin object
district_owin <- as.owin(boundary_district)
ppp_obj <- as.ppp(st_geometry(conflict_district))
# Handle duplicates
ppp_obj <- rjitter(ppp_obj, retry = TRUE, nsim = 1, drop = TRUE)
# Mask ppp object with the boundary
district_ppp_owin <- ppp_obj[district_owin]
district_ppp_owin <- rescale(district_ppp_owin, 1000, "km")
# Calculate the mean center
mean_coords <- c(mean(district_ppp_owin$x), mean(district_ppp_owin$y))
# Calculate the standard distance
sd_x <- sd(district_ppp_owin$x)
sd_y <- sd(district_ppp_owin$y)
standard_distance <- sqrt(sd_x^2 + sd_y^2)
# Prepare to plot for the current district
plot(district_ppp_owin,
main = paste("Standard Distance in", district),
pch = 19,
col = "lightgrey",
xlim = range(district_ppp_owin$x) + c(-5, 5), # Expand limits
ylim = range(district_ppp_owin$y) + c(-5, 5)) # Expand limits
# Highlight the mean center with an 'X'
points(mean_coords[1], mean_coords[2], pch = 4, col = "red", cex = 1) # Smaller 'X'
# Create a circle around the mean center for the standard distance
bearing <- seq(0, 2 * pi, length.out = 360)
circle_x <- mean_coords[1] + standard_distance * cos(bearing)
circle_y <- mean_coords[2] + standard_distance * sin(bearing)
# Draw the circle
lines(circle_x, circle_y, col = 'red', lwd = 2)
}
# Reset the plot layout
par(mfrow = c(1,1), mar = c(0,0,1,0)) 
5. 1st Order Spatial Point Patterns Analysis
5.1 Kernel Density Estimation
5.1.1 Working with Fixed Bandwidth Methods
Using the geospatial data sets prepared, I will now perform 1st order spatial point pattern analysis by leveraging kernel density estimation (KDE) to understand the intensity of conflicts in different regions.
I will be using a variety of fixed bandwidth methods via density() of the spatstat package, to determine the most optimal method for this analysis. Namely using bw.diggle(), bw.ppl(), bw.CvL() and bw.scott().
Steps taken to calculate the KDE:
- Extract the masked ppp object for the current quarter.
- Compute the kernel density estimate by setting the
signmaparameters.- Plot the kernel density estimate using
plot()where “Bw” represents the optimal bandwidth
For the purposes of identifying the most optimal bandwidth method, I will create a ppp_obj using the 2021 Q1 conflict events first to assist my decision-making.
# Set Up
ppp_obj = masked_ppp_list_km$`2021 Q1`
colours <- colorRampPalette(c("midnightblue", "skyblue"))(100)The bw.diggle() bandwidth is referred to as Diggle’s cross-validation bandwidth which minimises the mean-squared error (MSE) to balance between under and over-smoothing.
# bw.diggle()
kde_conflict_bw_diggle <- density(ppp_obj,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian")
optimal_bw_d = floor(bw.diggle(ppp_obj)[[1]]*10)/10
plot(kde_conflict_bw_diggle, main = paste("BW: diggle", "(",optimal_bw_d,"km)"), col = colours)
The second bandwidth method I attempted using is bw.ppl(), This method chooses the bandwidth that minimises the likelihood cross-validation score and improving the prediction accuracy of the kernel density estimate.
# bw.ppl()
kde_conflict_bw_ppl <- density(ppp_obj,
sigma=bw.ppl,
edge=TRUE,
kernel="gaussian")
optimal_bw_p = floor(bw.ppl(ppp_obj)[[1]]*10)/10
plot(kde_conflict_bw_ppl, main = paste("Bw: ppl", "(",optimal_bw_p,"km)"), col = colours)
Thirdly, let’s explore the bandwidth method bw.CvL(), also known as Cronie and Van Lieshout cross-validation, designed to provide an optimal, adaptive bandwidth for inhomogeneous point patterns.. Similar to bw.ppl(), it aims to reduce the error measure but also aims to balance over and under-fitting based on the spatial structure of the data.
# bw.CvL()
kde_conflict_bw_CvL <- density(ppp_obj,
sigma=bw.CvL,
edge=TRUE,
kernel="gaussian")
optimal_bw_c = floor(bw.CvL(ppp_obj)[[1]]*10)/10
plot(kde_conflict_bw_CvL, main = paste("Bw: CvL (",optimal_bw_c,"km)"), col = colours)
Lastly, I will explore the bw.scott() bandwidth method. This method returns separate bandwidths for the x- and y-axes which is ideal for our spatial data that contains both x and y components. I will combine these bandwidths into a single value for isotropic kernel density estimation by taking the taking the geometric mean as shown in the value returned by sigma_combined.
# bw.scott()
bw_values <- bw.scott(ppp_obj)
sigma_x <- bw_values[1]
sigma_y <- bw_values[2]
sigma_combined <- sqrt(sigma_x * sigma_y)
kde_conflict_bw_scott <- density(ppp_obj,
sigma = sigma_combined,
edge = TRUE,
kernel = "gaussian")
optimal_bw_s = floor(sigma_combined*10)/10
plot(kde_conflict_bw_scott, main = paste("Bw: scott", "(",optimal_bw_s,"km)"), col = colours)
Selecting a Bandwidth Method
Pramanik N. (2019) proposes that choosing an optimal bandwidth method is crucial in fitting the data appropriately to balance between bias and variance. As such, the four methods listed below cater to different types of data depending on how varied the densities are spread across and the granularity of conflict events. Let’s explore which bandwidth method is optimal for our dataset.
bw.diggle():appears to be effective for homogeneous data in seeing general conflict hotspots and uses the smallest bandwidth sizebw.ppl(): tends to choose slightly larger bandwidths, it provide more localised density estimates which highlights finer spatial details. As such, we can see more variability and finer details in the density distribution, with more variation between high- and low-density areas.bw.CvL(): the kernel density plot shows that CvL makes a good attempt in balancing between detail and smoothness, making it more suitable for capturing the overall density trends in spatial data with some local structures highlighted.bw.scott():as shown, the geometric mean ensures equal smoothing in both x and y directions, and is largely similar tobw.Cvl(), making it a good choice for a balanced and general and fast overview of the spatial data distribution.
Plot all bandwidth methods
par(mfrow = c(2,2), mar = c(0,0,1,0))
plot(kde_conflict_bw_diggle, main = paste("diggle (",optimal_bw_d,"km)"), col = colours)
plot(kde_conflict_bw_ppl, main = paste("ppl (",optimal_bw_p,"km)"), col = colours)
plot(kde_conflict_bw_CvL, main = paste("CvL (",optimal_bw_c,"km)"), col = colours)
plot(kde_conflict_bw_scott, main = paste("scott (",optimal_bw_s,"km)"), col = colours)
# Plot histogram to compare KDE
par(mar = c(2,2,2,2),mfrow = c(2,2))
hist(kde_conflict_bw_diggle, main = "diggle")
hist(kde_conflict_bw_ppl, main = "ppl")
hist(kde_conflict_bw_CvL, main = "scott")
hist(kde_conflict_bw_scott, main = "CvL")
💡 Decision: I decided to use
bw.CvL()for computing the KDE of the masked ppp objects based on each quarter. The KDE function outputs a relatively smooth density estimate that isn’t too detailed likebw.ppl()and not as generalised asbw.scott().This coincides with wider range of spatial concentrations captured by the histogram, doing a good job with both capturing the non-homogeneous data spread of our Myanmar conflict data. Moreover, it isn’t as computationally heavy as
bw.ppl().
Putting Together our Fixed KDE using bw.CvL()
Now, let us perform the KDE computation for the conflict events across all quarters using bw.CvL().
# Calculate density using bw.CvL()
par(mfrow = c(2,3), mar = c(0,0,1,0))
for (quarter in names(masked_ppp_list_km)) {
ppp_obj = masked_ppp_list_km[[quarter]]
kde_conflict_bw <- density(ppp_obj,
sigma=bw.CvL,
edge=TRUE,
kernel="gaussian")
optimal_bw = floor(bw.CvL(ppp_obj)[[1]]*10)/10
plot(kde_conflict_bw, main = paste(quarter, "(Bw:",optimal_bw,"m)"), col = colours)
}
Observations
Bandwidth Size: A bandwidth of around 60,000 to 100,000 is considered relatively large as compared to the bandwidth returned from using bw.diggle() and bw.ppl(). Hence, this results in a smoother density estimate with less emphasis on local clusters as indicated in the generalised spatial trends.
Note: Using a fixed KDE is beneficial for our quarterly analysis but the bandwidth returned is different for each quarter. Later on, I’ll take the average bandwidth across different time periods so we can focus on spatial distribution over time.
5.1.2 Working with Different Kernel Methods
Bias is an inevitable issue when it comes to choosing the type of kernel function. Rajagopalan B., et al (1997) states that kernels such as Epanechnikov, Quartic and Disc suffer a “leakage problem” where some of the probability mass is truncated at the boundary, leading to boundary bias. At the same time, the Gaussian kernel has no truncation, but small “leakage” beyond data range of the boundary inadvertantly causes bias to exist still. Shen et al. (2020) further proposes that the bandwidth is a more of an influential factor in kernel estimation than the selection of different kernel functions.
Let us evaluate if this is true by experimenting with a variety of kernel methods as they ultimately control how we weight points within the bandwidth radius. The default kernel in density.ppp() is the gaussian. alternatives such as epanechnikov, quartic, and disc are also available. I will use the 2021 Q1 conflict data to assist in identifying the most optimal kernel method.
# Set Up
par(mfrow = c(2,2), mar = c(0,0,1,0))
ppp_obj = masked_ppp_list_km$`2021 Q1`
# Using the gaussian kernel
kde_conflict_g <- density(ppp_obj,
sigma=bw.CvL,
edge=TRUE,
kernel="gaussian")
plot(kde_conflict_g, main="Gaussian Kernel", col = colours)
# Using the epanechniko kernel
kde_conflict_e <- density(ppp_obj,
sigma=bw.CvL,
edge=TRUE,
kernel="epanechnikov")
plot(kde_conflict_e, main="Epanechnikov Kernel", col = colours)
# Using the quartic kernel
kde_conflict_q <- density(ppp_obj,
sigma=bw.CvL,
edge=TRUE,
kernel="quartic")
plot(kde_conflict_q, main="Quartic Kernel", col = colours)
# Using the disc kernel
kde_conflict_d <- density(ppp_obj,
sigma=bw.CvL,
edge=TRUE,
kernel="disc")
plot(kde_conflict_d, main="Disc Kernel", col = colours)
# Plot histogram to compare KDE
par(mar = c(2,2,2,2),mfrow = c(2,2))
hist(kde_conflict_g, main = "Gaussian")
hist(kde_conflict_e, main = "Epanechnikov")
hist(kde_conflict_q, main = "Quartic")
hist(kde_conflict_d, main = "Disc")
Observations
As expected, we don’t see major variations in the smoothness and spread of KDE across a range of distances, and bias continues to be present across all kernel methods as seen from the long-tails. However, there are slight differences in the sharpness of the bandwidth radius in terms of how localised or widespread our data points are being captured.
- Gaussian: provides a localised density estimate over the entire spatial extent as compared to
epanechnikovandquartic. It is good at highlighting variance and opposing ends of conflict intensities as shown by the wider range used in the legend. - Epanechnikov: It is more efficient than the
gaussianin terms of variance but produces a slightly rougher surface. It is also more localised than thequartickernel, focusing on areas near each point, with a sharper boundary at the bandwidth limit. - Quartic: Results in a good balance between smoothness and localised influence, smoother than
epanechnikovbut with similar properties. It appears suitable for moderate smoothing and sharper focus on local patterns. - Disc: results in the sharpest density estimate as compared to the other three kernels as all points within a certain distance are made to have equal influence and zero influence beyond that distance.
Decision: Hence, I will use the quartic kernel method to ensure a relatively smooth density estimate with emphasis on local points over distant ones.
5.1.3.1 Using kernel = ‘quartic’
As such, I run the density estimate computation using kernel = 'quartic'.
# Using 'quartic' kernel
par(mfrow = c(2,3), mar = c(0,0,1,0))
for (quarter in names(masked_ppp_list_km)) {
ppp_obj = masked_ppp_list_km[[quarter]]
kde_conflict_bw <- density(ppp_obj,
sigma=bw.CvL,
edge=TRUE,
kernel="quartic")
optimal_bw = floor(bw.CvL(ppp_obj)[[1]]*1000)/1000
plot(kde_conflict_bw, main = paste(quarter, "(Bw:",optimal_bw,"km)"), col = colours)
}
Observations
We can see high densities of armed conflict in the central and southern regions of Myanmar but more can be uncovered from conflict data. Let’s proceed to the next section.
5.1.3.2 Using sigma = 71.831
For all subsequent fixed KDE computations, I will assign sigma using the average of the CvL bandwidth returned from each quarter. Here’s the calculations based on the plots returned:
Average bandwidth size =
61.649 + 64.386 + 74.501 + 114.08 + 103.863 + 103.863 + 95.567 + 56.757 + 57.752 + 66.968 + 49.649 + 48.135 + 48.135 + 60.323) / 14 = 71.831
Let’s recompute the Fixed KDE based on the newly calculated average bandwidth such that sigma = 71.831. I’ll store the quarterly KDE outputs into a list called kde_conflict_bw_list.
# Add KDE into this list
kde_conflict_bw_list <- list()
for (quarter in names(masked_ppp_list_km)) {
ppp_obj = masked_ppp_list_km[[quarter]]
kde_conflict_bw <- density(ppp_obj,
sigma=71.831,
edge=TRUE,
kernel="quartic")
kde_conflict_bw_list[[quarter]] <- kde_conflict_bw
}# Plot graph
par(mfrow = c(2,3), mar = c(0,0,1,0))
for (quarter in names(kde_conflict_bw_list)){
kde_conflict_bw <- kde_conflict_bw_list[quarter]
plot(kde_conflict_bw, main = paste(quarter, "(Bw: 71.831 km)"), col = colours)
}
5.1.3 Working with Adaptive KDE
As seen above, fixed bandwidths tend to oversmooth the mode of the distribution. On the contrary, the adaptive kernel estimate has the ability to reduce variability of estimates in areas with low density and increases it in areas with higher density (The Stata Journal, 2003).
Once again, let us use the 2021 Q1 conflict data to illustrate the difference in outputs of all three adaptive methods.
# Set Up
ppp_obj = masked_ppp_list_km$`2021 Q1`
# Using Voronoi adaptive KDE
vd_adaptive_kde <- adaptive.density(ppp_obj, f=1, method="voronoi")
# Plot
par(mar = c(0,1,1,1))
plot(vd_adaptive_kde, main = "Voronoi-Dirichlet Adaptive KDE", col = colours)
# Using adaptive KDE
adaptive_kde <- adaptive.density(ppp_obj, method="kernel")
# Plot
par(mar = c(0,1,1,1))
plot(adaptive_kde, main = "Adaptive KDE", col = colours)
I’ve used a relatively larger number of neighbours (i.e. k = 10) to provide a smoother, more general density estimate to capture broader trends and may smooth out details.
# Using nearest neighbour adaptive KDE
nn_kde <- nndensity(ppp_obj, k=10)
# Plot
par(mar = c(0,1,1,1))
plot(nn_kde, main = "Nearest Neighbour Adaptive KDE", col = colours)
We can also compare the performance of each method based on the top 4 states with highest proportion of conflicts as highlighted earlier.
Comparing the three Adaptive KDE Types
From the outputs above, it appears that there is no major differences between the distribution of KDE values returned across the three methods, where there is high concentration of points in a specific area. Hence, we will choose to go with Adapative Kernel method.
Let’s compare the results of my two selected fixed and adaptive KDEs (E.g. Magway District)
We can observe how adaptive kernels provides a more detailed picture of conflict spatial distribution but since it’s largely localised, conflict spots require more effort in identifying and can be computationally heavy for this exercise.
Additionally, varying bandwidth makes comparisons across regions or time periods (like quarters) more difficult because the scale of smoothing is not constant across space and time.
5.1.4 Converting Gridded KDE Output into Raster
Next, we need to convert the KDE output to KDE raster layers before it can be viewed using tmap.
Step 1) Converting KDE to Spatial Grid Data Frame
library(grid)
plot_list <- list()
for (quarter in names(kde_conflict_bw_list)) {
ppp_obj <- kde_conflict_bw_list[[quarter]]
gridded_ppp_obj <- as(ppp_obj, "SpatialGridDataFrame")
plot_list[[quarter]] <- spplot(gridded_ppp_obj, main = paste(quarter), col.regions = colours)
}
library(gridExtra)
plot_list_subset1 <- plot_list[1:6]
plot_list_subset2 <- plot_list[7:12]
plot_list_subset3 <- plot_list[13:14]
grid.newpage()
grid.arrange(grobs = plot_list_subset1, ncol = 3, nrow = 2)
grid.newpage()
grid.arrange(grobs = plot_list_subset2, ncol = 3, nrow = 2)
grid.newpage()
grid.arrange(grobs = plot_list_subset3, ncol = 3, nrow = 2)
Step 2) Rasterisation of Grid Outputs & Assigning Projection Systems
gridded_ppp_obj_raster_list <- list()
for (quarter in names(kde_conflict_bw_list)) {
gridded_ppp_obj = kde_conflict_bw_list[[quarter]]
gridded_ppp_obj_raster <- raster(gridded_ppp_obj)
projection(gridded_ppp_obj_raster) <- CRS("+init=EPSG:32647")
gridded_ppp_obj_raster_list[[quarter]] <- gridded_ppp_obj_raster
}
# Inspect for 2024 Q2
gridded_ppp_obj_raster_list$`2024 Q2`class : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 7.302001, 16.30032 (x, y)
extent : -210.0086, 724.6476, 1072.026, 3158.467 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=47 +datum=WGS84 +units=m +no_defs
source : memory
names : layer
values : 7.061861e-06, 0.01847887 (min, max)Step 3) Plot Maps
tmap_mode("plot")
plots_by_quarter <- list()
for (quarter in names(gridded_ppp_obj_raster_list)){
gridded_ppp_obj_raster = gridded_ppp_obj_raster_list[[quarter]]
raster_plot <- tm_shape(gridded_ppp_obj_raster) +
tm_raster("layer", title="Density", palette="Blues") +
tm_layout(legend.position = c("left","bottom"),frame=FALSE, main.title = quarter,
main.title.size=1, main.title.position = "center", legend.text.size = 0.5,legend.title.size = 0.7)
plots_by_quarter[[quarter]] <- raster_plot
}
tmap_arrange(plots_by_quarter[1:14], ncol=5, nrow=3)
Observations
Plotting raster grid versions of KDE outputs uses discrete colour ranges which does a good job in highlighting gradual changes in conflict events across an area. Since 2021 Q2, more conflicts are seen in Southern parts of Myanmar. The density range differs for each quarter but we can see an increase in no. of armed conflicts per kilometre from 2021 Q1 to 2022 Q2, which stagnates in density and increases again in 2023 Q4.
5.2 Nearest Neighbour Analysis
Our current analyses does not reveal patterns of clustering or dispersion, to which Michael J. Crawley proposes to employ Clark-Evans test spatial randomness for its simplicity and applicability for first-order spatial analysis, which means checking for overall spatial randomness based on nearest-neighbor distances. (Crawley M. J. , 2007)
Clark-Evans Test
The test checks whether the observed point pattern of armed conflicts in Myanmar shows clustering (points are closer than expected under randomness), dispersion (points are more spread out), or randomness.
5.2.1 Clark-Evans Test (Myanmar)
The test hypotheses are:
- Ho = The distribution of armed conflicts in Myanmar are randomly distributed.
- H1= The distribution of armed conflicts in Myanmar are not randomly distributed.
- The 95% confident interval will be used.
We will conduct the test using clarkevans.test() of statspat.
for (quarter in names(masked_ppp_list_km)) {
ppp_obj = masked_ppp_list_km[[quarter]]
print(quarter)
print(clarkevans.test(ppp_obj,
correction="none",
clipregion="boundary_sf",
alternative=c("clustered"),
nsim=99))
}[1] "2024 Q2"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.26663, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2024 Q1"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.23563, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2023 Q4"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.21795, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2023 Q3"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.22002, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2023 Q2"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.24485, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2023 Q1"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.24365, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2022 Q4"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.22139, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2022 Q3"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.23974, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2022 Q2"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.22989, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2022 Q1"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.21976, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2021 Q4"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.21341, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2021 Q3"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.21808, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2021 Q2"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.17458, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
[1] "2021 Q1"
Clark-Evans test
No edge correction
Z-test
data: ppp_obj
R = 0.24696, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
Observations
For a 95% confidence level, If the p-value < 0.05, I will reject the null hypothesis of complete spatial randomness and check if data is uniform (R > 1) or clustered (R < 1).
With that said, all tests conducted across each quarter rejects the null hypothesis as p < 0.05 and spatial points are found to be clustered since R < 1.
5.2.2 Clark-Evans Test (Top 4 Districts)
# Yinmarbin District
clarkevans.test(yinmarbin_ppp_owin,
correction="none",
clipregion="boundary_yinmarbin",
alternative=c("clustered"),
nsim=39)
Clark-Evans test
No edge correction
Z-test
data: yinmarbin_ppp_owin
R = 0.13718, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)# Shwebo District
clarkevans.test(shwebo_ppp_owin,
correction="none",
clipregion="boundary_shwebo",
alternative=c("clustered"),
nsim=39)
Clark-Evans test
No edge correction
Z-test
data: shwebo_ppp_owin
R = 0.20141, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)# Pakokku District
clarkevans.test(pakokku_ppp_owin,
correction="none",
clipregion="boundary_pakokku",
alternative=c("clustered"),
nsim=39)
Clark-Evans test
No edge correction
Z-test
data: pakokku_ppp_owin
R = 0.20698, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)# Mandalay District
clarkevans.test(mandalay_ppp_owin,
correction="none",
clipregion="boundary_mandalay",
alternative=c("clustered"),
nsim=39)
Clark-Evans test
No edge correction
Z-test
data: mandalay_ppp_owin
R = 0.16014, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
Observations
I will reject the null hypothesis of complete spatial randomness since the p values of each district is smaller than 0.05. Additionally, we can observe clustering of spatial points since R < 1 is returned for all districts. However, this test isn’t sufficient in highlighting the statistical significance of the spatial patterns - I’ll use the Monte Carlo Simulation later to handle this.
5.3 Further Data Exploration
By using the fixed KDE with CvL bandwidth and quartic kernel, let’s see what insights can we glean from the density of conflicts in Myanmar.
5.3.1 KDE by Event Type
First, let’s identify the unique event types in this dataset.
# Check unique events
unique(conflict_data_sf$event_type)[1] "Battles" "Strategic developments"
[3] "Violence against civilians" "Explosions/Remote violence"Now, let us analyse the kernel density estimate of each unique event type found in conflict_data_sf to identify hot and cold spots across Myanmar.
Plot the KDE based on Event Type
# Set Up
par(mfrow = c(2,2), mar = c(0,0,1,0))
conflict_data_sf %>%
group_by(event_type) %>%
group_split() -> conflict_by_event_type
# Convert the sf object to owin
district_boundary <- as.owin(st_as_sfc(boundary_sf))
kde_list <- lapply(seq_along(conflict_by_event_type), function(i) {
data <- conflict_by_event_type[[i]]
event_type <- unique(data$event_type)
ppp_obj <- as.ppp(st_geometry(data), W = district_boundary)
ppp_obj <- rescale(ppp_obj, 1000, "km")
kde <- density(ppp_obj,
sigma=71.831,
edge=TRUE,
kernel="quartic")
plot(kde, main = paste(event_type), col=colours)
return(kde)
})
Observations
We can almost see an equal spread of all four event types, with explosions and strategic violence being more dominantly found in Central-Western Myanmar, in the Sagaing and Mandalay states, followed by battles and violence against citizens. All events seem to plaque the Southern states (e.g. Yangon) with the exception of battles.
5.3.2 KDE Across Top 4 States With Most Conflicts
Previously, we identified the top 4 states with the highest proportions of conflicts as Sagaing, Mandalay, Magway and Yangon. We can delve deeper into each state by analysing the intensity of conflicts across these states using density().
Plot the KDE of Top 4 States
# Set Up
par(mfrow = c(2,2), mar = c(0,0,1,0))
# Sagaing
district_boundary <- as.owin(st_as_sfc(boundary_sagaing))
ppp_obj <- as.ppp(st_geometry(conflict_sagaing), W = district_boundary)
ppp_obj_sagaing <- rescale(ppp_obj, 1000, "km")
kde <- density(ppp_obj_sagaing,
sigma=71.831,
edge=TRUE,
kernel="quartic")
plot(kde, main = paste("Sagaing"), col=colours)
# Mandalay
district_boundary <- as.owin(st_as_sfc(boundary_mandalay))
ppp_obj <- as.ppp(st_geometry(conflict_mandalay), W = district_boundary)
ppp_obj_mandalay <- rescale(ppp_obj, 1000, "km")
kde <- density(ppp_obj_mandalay,
sigma=71.831,
edge=TRUE,
kernel="quartic")
plot(kde, main = paste("Mandalay"), col=colours)
# Magway
district_boundary <- as.owin(st_as_sfc(boundary_magway))
ppp_obj <- as.ppp(st_geometry(conflict_magway), W = district_boundary)
ppp_obj_magway <- rescale(ppp_obj, 1000, "km")
kde <- density(ppp_obj_magway,
sigma=71.831,
edge=TRUE,
kernel="quartic")
plot(kde, main = paste("Magway"), col=colours)
# Yangon
district_boundary <- as.owin(st_as_sfc(boundary_yangon))
ppp_obj <- as.ppp(st_geometry(conflict_yangon), W = district_boundary)
ppp_obj_yangon <- rescale(ppp_obj, 1000, "km")
kde <- density(ppp_obj_yangon,
sigma=71.831,
edge=TRUE,
kernel="quartic")
plot(kde, main = paste("Yangon"), col=colours)
Observations
It’s interesting that armed conflict isn’t evenly distributed across the states though it does seem that armed conflict has inflicted the entire state of Yangon. Nonetheless, it is worth noting that Yangon is relatively smaller in size than the other three states and that will increase the density of conflict quite significantly.
5.3.3 KDE of Top 4 States by Event Type
It’ll also be interesting to breakdown each top 4 state by the event type category as shown.
Observations
The kernel density of violence against civillians is generally found to be the lowest amongst all conflict events. Additionally, all types of armed conflicts tend to occur repeatedly in the same parts of each state. E.g. conflicts regarding strategic development tend to happen in Southern part of the Sagaing state, just as it is for explosions/remote violence.
5.5.4 KDE by Interaction Type
# Convert Interaction Code to Text
unique(conflict_data_sf$interaction) [1] 13 10 12 37 70 17 22 11 33 30 27 20 23 60 47 28 38 80 78 18 24 14library(dplyr)
# Create a named vector for mapping interaction values
interaction_map <- c(
"10" = "SOLE STATE FORCES ACTION",
"11" = "STATE FORCES VERSUS STATE FORCES",
"12" = "STATE FORCES VERSUS REBELS",
"13" = "STATE FORCES VERSUS POLITICAL MILITIA",
"14" = "STATE FORCES VERSUS IDENTITY MILITIA",
"17" = "STATE FORCES VERSUS CIVILIANS",
"18" = "STATE FORCES VERSUS EXTERNAL/OTHER FORCES",
"20" = "SOLE REBEL ACTION",
"22" = "REBELS VERSUS REBELS",
"23" = "REBELS VERSUS POLITICAL MILITIA",
"24" = "REBELS VERSUS IDENTITY MILITIA",
"27" = "REBELS VERSUS CIVILIANS",
"28" = "REBELS VERSUS OTHERS",
"30" = "SOLE POLITICAL MILITIA ACTION",
"33" = "POLITICAL MILITIA VERSUS POLITICAL MILITIA",
"34" = "POLITICAL MILITIA VERSUS IDENTITY MILITIA",
"37" = "POLITICAL MILITIA VERSUS CIVILIANS",
"38" = "POLITICAL MILITIA VERSUS OTHERS",
"40" = "SOLE IDENTITY MILITIA ACTION",
"44" = "IDENTITY MILITIA VERSUS IDENTITY MILITIA",
"47" = "IDENTITY MILITIA VERSUS CIVILIANS",
"48" = "IDENTITY MILITIA VERSUS OTHER",
"70" = "SOLE CIVILIAN ACTION",
"77" = "CIVILIANS VERSUS CIVILIANS",
"78" = "OTHER ACTOR VERSUS CIVILIANS",
"80" = "SOLE OTHER ACTION",
"88" = "OTHER VERSUS OTHER"
)
conflict_data_sf <- conflict_data_sf %>%
rename(interaction_code = interaction) %>%
mutate(interaction = interaction_map[as.character(interaction_code)])
# View the updated dataframe
head(conflict_data_sf[c("interaction_code","interaction")])Simple feature collection with 6 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 76997.72 ymin: 2428487 xmax: 214961 ymax: 2533434
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 6 × 3
interaction_code interaction geometry
<dbl> <chr> <POINT [m]>
1 13 STATE FORCES VERSUS POLITICAL MILI… (214961 2452068)
2 13 STATE FORCES VERSUS POLITICAL MILI… (198303.2 2499463)
3 13 STATE FORCES VERSUS POLITICAL MILI… (189105.4 2533434)
4 10 SOLE STATE FORCES ACTION (160913.9 2522331)
5 13 STATE FORCES VERSUS POLITICAL MILI… (146213 2428487)
6 10 SOLE STATE FORCES ACTION (76997.72 2447719)1. Involvement of Political Militia
The conflict landscape in Myanmar is complex and fragmented, with political militias focused on advancing political agendas or controlling governance structures. Let’s discover how the density of armed conflicts differ for each interaction of the political militia. They can be associated as either of these groups.
- Pro-Government Militias: supporting the central government or Tatmadaw, such as Pyithu Sit and BGFs.
- Pro-Democracy Militias: emerged as part of the resistance movement against the military junta. E.g. the People’s Defence Forces (PDF)
Plot conflicts involving political militia
codes <- c('13','23','33')
par(mfrow = c(3,1), mar = c(0,0,1,0))
for(code in codes) {
# Filter for the specific interaction type
filtered_data <- conflict_data_sf %>%
filter(interaction_code == code)
# Convert to ppp object
boundary_owin <- as.owin(boundary_sf)
ppp_owin_obj <- as.ppp(filtered_data$geometry, W = boundary_owin)
ppp_owin_obj_km <- rescale(ppp_owin_obj, 1000, "km")
# Compute KDE
kde <- density(ppp_owin_obj_km,
sigma=71.831,
edge=TRUE,
kernel="quartic")
# Plot the KDE
plot(kde, main = paste(interaction_map[code]), col = colours)
}
Observations
The political army or organisations tend to be less involved with unarmed civilians but state forces, rebels and other policial militia.
- Against State Forces
- High density of conflicts with state forces (e.g. police) in central and western Myanmar
- These regions have significant ethnic diversity and political grievances. E.g. Rakhine State in the west is home to ethnic Rakhine groups who have historically sought greater autonomy and have been involved in conflicts with state forces.
- Central regions, such as Sagaing and Magway, have also seen intense resistance against state forces, particularly after the 2021 military coup.
- Against Rebels
- Rebels and political militias operate throughout Myanmar with varying objectives but of a much smaller intensity than with other actors (seen from its low density scale).
- Clashes between rebels and political militias can stem from differing political goals, territorial disputes, or strategic interests.
- For instance, a rebel group might view a political militia as a rival in a contested region or as a competitor for local support.
- Against other Political Miliitia
- In central and western Myanmar, where political militias are more active, there is often competition among different resistance groups. These militias may have different leadership, ideologies, or goals, leading to clashes.
2. Involvement of Identity Militia
Identity militias are typically organised around ethnic, religious, or cultural identities. heir main goal is to protect the interests, rights, and autonomy of specific identity-based groups within Myanmar.
Plot conflicts involving identity militia
codes <- c('14','24')
par(mfrow = c(2,1), mar = c(0,0,1,0))
for(code in codes) {
# Filter for the specific interaction type
filtered_data <- conflict_data_sf %>%
filter(interaction_code == code)
# Convert to ppp object
boundary_owin <- as.owin(boundary_sf)
ppp_owin_obj <- as.ppp(filtered_data$geometry, W = boundary_owin)
ppp_owin_obj_km <- rescale(ppp_owin_obj, 1000, "km")
# Compute KDE
kde <- density(ppp_owin_obj_km,
sigma=71.831,
edge=TRUE,
kernel="quartic")
# Plot the KDE
plot(kde, main = paste(interaction_map[code]), col = colours)
}
Observations
Armed conflicts involving the identity militia is of a much smaller intensity than political militia with about 5e^-0.5 or 3 conflict events per kilometre.
- State Forces vs Identity Militia
Conflicts between these two parties are concerntrated in the Western region which has a history of conflict, including uprisings against colonial and post-colonial governments.
Additionally, the region’s proximity to Bangladesh has some influence on local conflicts as ethnic groups in Myanmar may have connections or support from across the border.
In the Rakhine State of western Myanmar, historical grievances and demands for autonomy has also fueled tensions between these groups and the central government (Crisis Group, 2024).
- Rebels vs identity militia
- Higher intensity of armed conflicts are seen in North Myanmar
- Northern Myanmar hosts various ethnic armed groups and rebel factions that seek autonomy or independence. (E.g. Kachin, Shan and Chin)
- Northern Myanmar is rich in natural resources, e.g. jade and timber. Control over these resources can be a significant driver of conflict, as various armed groups vie for control and when there’s a lack of economic opportunities. (Fishbein & Lusan, 2022)
- Unresolved conflicts from colonial and post-colonial periods continue to influence current tensions between different ethnic groups and armed factions in the Northen region.
3. Involvement of Civilian Actors
Plot conflicts involving civilians
codes <- c('17', '27', '37', '47', '78')
par(mfrow = c(3,2), mar = c(0,0,1,0))
for(code in codes) {
# Filter for the specific interaction type
filtered_data <- conflict_data_sf %>%
filter(interaction_code == code)
# Convert to ppp object with boundary window
boundary_owin <- as.owin(boundary_sf)
ppp_obj <- as.ppp(filtered_data$geometry, W = boundary_owin)
# Rescale the ppp object
ppp_obj_km <- rescale(ppp_obj, 1000, "km")
# Compute KDE
kde <- density(ppp_obj_km,
sigma = 71.831,
edge = TRUE,
kernel = "quartic")
# Plot the KDE
plot(kde, main = paste(interaction_map[code]), col = colours)
}
Observations
Across all actors, we see that perhaps civilians are involved with the biggest range of actors.
Notably, we see the lowest intensity of clashes with identity militia which could stem from how civilians might be part of the same social/ethnic networks as the identity militia and these militia may seek the support of civilian populations too.
In contrast, state forces and civilians have the highest intensity of armed conflicts
In Western Rakhine State, Buddhist-majority and Muslim-minority communities have experienced significant tension causing state forces to act in response to these tensions which has led to severe clashes with civilian populations.
Reports of human rights abuses, including arbitrary arrests, torture, and extrajudicial killings by state forces, has led to large-scale displacement and humanitarian crises, further exacerbating tensions between state forces and civilians. The struggle for resources and safety often intensifies the conflicts. (Farge E. & Mantovani C., 2024)
4. Civilians Involvement by Event Type
Compute KDE of Civilians by Event Type
# Set Up
par(mfrow = c(2,2), mar = c(0,0,1,0))
codes <- c('17', '27', '37', '47', '78')
conflict_data_sf %>%
filter(interaction_code %in% codes) %>%
group_by(event_type) %>%
group_split() -> conflict_by_event_type
# Convert the sf object to owin
district_boundary <- as.owin(st_as_sfc(boundary_sf))
kde_list <- lapply(seq_along(conflict_by_event_type), function(i) {
data <- conflict_by_event_type[[i]]
event_type <- unique(data$event_type)
ppp_obj <- as.ppp(st_geometry(data), W = district_boundary)
ppp_obj <- rescale(ppp_obj, 1000, "km")
kde <- density(ppp_obj,
sigma=71.831,
edge=TRUE,
kernel="quartic")
plot(kde, main = paste(event_type), col=colours)
return(kde)
})
Observations
Civilians are seen to be embroiled mostly in conflicts resulting from strategic development and violence against civilians events, primarily in Central and Western Myanmar. Fortunately, conflicts involving explosions or remote violence is not as intense against civilians but this raises a great concern on the humanitarian crisis faced by civilians in Myanmar.
5.5.5 OpenStreetMap Myanmar - Spatial Points
Using tmap functions, I will display an interactive view of te KDE layers on openstreetmap of Myanmar to observe intensity of conflicts at the district level. We can see higher intensities of conflict in central districts (e.g. Yinmarbin, Pakokku and Shwebo) and western districts (e.g. Yangon).
ppp_obj<- as.ppp(conflict_data_sf$geometry)
ppp_owin_obj <- ppp_obj[myanmar_owin]
ppp_owin_obj_km <- rescale(ppp_owin_obj, 1000, "km")
kde_fixed <- density(ppp_owin_obj_km, sigma=bw.CvL, edge=TRUE, kernel="quartic")
raster_kde_fixed <- raster(kde_fixed)
projection(raster_kde_fixed) <- CRS("+init=EPSG:32647 +units=km")# Plot KDE Map on OpenStreetMap
tmap_mode('view')tmap mode set to interactive viewingkde_fixed_output <- tm_basemap(server = "OpenStreetMap.HOT") +
tm_basemap(server = "Esri.WorldImagery") +
tm_shape(raster_kde_fixed) +
tm_raster("layer",
n = 10,
title = "KDE Fixed CvL",
alpha = 0.6,
palette = c("#fafac3","#fd953b","#f02a75","#b62385","#021c9e")) +
tm_shape(boundary_sf)+
tm_polygons(alpha=0.1,id="DT")+
tmap_options(check.and.fix = TRUE)
kde_fixed_outputtmap_mode("plot")tmap mode set to plotting6. 2nd Order Spatial Point Patterns Analysis
Unlike 1st-order analysis, which studies the intensity of points (e.g., density), let’s also leverage 2nd-order analysis to examine how points are distributed relative to each other, which can offer deeper insights into the spatial interaction between events.
I will use the K-function and L-function is to understand the spatial relationships between events, particularly focusing on whether the points exhibit clustering, uniformity, or randomness.
6.1 Using K-Function Estimation
K-function helps detect spatial patterns by comparing the observed distribution of points against a random pattern at different distances.
6.1.1 Yinmarbin District
1) Computing K-Function Estimation
For Yinmarbin district, let’s compute K-function estimates by using Kest() of the spatstat package by using the yinmarbin_ppp_owin masked ppp object we created previously.
K_ck = Kest(yinmarbin_ppp_owin, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(km)",
main = paste("Yinmarbin District (K-Function)")).png)
Observations
How to interpret the plot:
- K-iso represents the observed or estimated K-function value calculated from the actual data
- K-pois is the theoretical K-function that represents the expected K-function
With that said…
We can observe how the observed line (K-iso) is found to lie above the theoretical line (K-pois) which suggests conflict points in Yinmarbin are highly clustered. Clustering is the strongest at a 20km distance which suggests large-scale clustering. In fact, it is more clustered together than expected by the null hypothesis.
Note: Since I had used the default edge = TRUE settings, edge correction will account for missing neighbours outside the boundary which helps maintain an accurate estimate of the K-function. Hence, there is marginal difference in the actual and expected K-function.
2) Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test (i.e. Monte Carlo simulation test) will be conducted. The hypothesis and test are as follows:
- Ho = The distribution of conflict events in Myanmar are randomly distributed.
- H1= The distribution of conflict events in Myanmar are not randomly distributed.
- The null hypothesis will be rejected if the observed K-function lies above/below the theoretical K-function and envelope.
By using envelope(), we can get a more robust interpretation by comparing the observed K-function against a simulation envelope of K-functions generated under the null hypothesis.
Note
To achieve a 95% confidence envelope in a K-function test with Complete Spatial Randomness, I will need to exclude the upper 2.5% and lower 2.5% of the simulated K-functions., i.e. I will need to generate at least 40 simulations where nsim = 39.
We’ll set rank = 1 for a conservative setting, such that the envelope is based on the extreme values (highest and lowest) from the simulations.
# Monte Carlo test with K-function
K_yinmarbin_csr <- envelope(yinmarbin_ppp_owin, Kest,
nsim = 39, rank = 1, glocal=TRUE)
plot(K_yinmarbin_csr, main = paste("Yinmarbin District (CSR)"))
Observations
The tight envelope suggests we can be confident that deviations between the observed and theoretical lines are meaningful, rather than due to random variation. It also indicates that simulations under CSR are producing homogeneous patterns.
On the other hand, the observed K-function line is constantly above the shaded region of the envelope and theoretical line, suggesting that observed spatial pattern are more clustered than expected under the null hypothesis of CSR
Just to be sure that there is no improvement in a higher number of simulations, I ran it again with nsim = 99. As shown below, the k-function output is similar but with a somewhat smoother plot than before since we are working with less variability. Hence, I’ll stick to nsim = 39.
6.1.2 Shwebo District
1) Computing K-Function Estimation
For Shwebo district, let’s compute K-function estimates by using Kest() of the spatstat package.
K_ck = Kest(shwebo_ppp_owin, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(km)",
main = paste("Shwebo District (K-Function)"))The observed line is consistently above at all distances, implying that conflict events are clustered at both small and large areas, rather than being randomly distributed across the area.
 (1).png)
Observations
The observed line is consistently above at all distances, implying that conflict events are clustered at both small and large areas, rather than being randomly distributed across the area.
2) Performing Complete Spatial Randomness Test
# Monte Carlo test with K-function
K_shwebo_csr <- envelope(shwebo_ppp_owin, Kest,
nsim = 39, rank = 1, glocal=TRUE)
plot(K_shwebo_csr, main = paste("Shwebo District (CSR)"))
Observations
Since the observed line lies above the envelope and the shaded region of the envelope is rather narrow, it strongly suggests clustering in Shwebo, rather than due to random variation.
6.1.3 Pakokku District
1) Computing K-Function Estimation
For Pakokku district, let’s compute K-function estimates by using Kest() of the spatstat package.
K_ck = Kest(pakokku_ppp_owin, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(km)",
main = paste("Pakokku District (K-Function)"))
Observations
We can observe how the observed line (K-iso) is found to be above the theoretical line (K-pois) which suggests conflict points in Pakokku are highly clustered. In fact, it is more clustered together than expected by the null hypothesis.
2) Performing Complete Spatial Randomness Test
# Monte Carlo test with K-function
K_pakokku_csr <- envelope(pakokku_ppp_owin, Kest,
nsim = 39, rank = 1, glocal=TRUE)
plot(K_pakokku_csr, main = paste("Pakokku District (CSR)"))
Observations
Likewise, we can be certain about the clustering patterns in Pakokku since the observed K values deviates above the envelope as shown.
6.1.4 Mandalay District
1) Computing K-Function Estimation
For Mandalay district, let’s compute K-function estimates by using Kest() of the spatstat package.
K_ck = Kest(mandalay_ppp_owin, correction = "Ripley")
plot(K_ck, . -r ~ r, ylab= "K(d)-r", xlab = "d(km)",
main = paste("Mandalay District (K-Function)"))
Observations
Significant clustering is observed generally across multiple distances in Mandalay and clusterinng is particularly evident at longer spatial distances in Mandalay seen from how the deviations of observed K values increase by the distance.
2) Performing Complete Spatial Randomness Test
# Monte Carlo test with K-function
K_mandalay_csr <- envelope(mandalay_ppp_owin, Kest,
nsim = 39, rank = 1, glocal=TRUE)
plot(K_mandalay_csr, main = paste("Mandalay District (CSR)"))
Observations
The envelope size continues to be narrow with K values plotted above the envelope. This strongly suggests that the point pattern is not random, and the points are clustered.
6.2 Using L-Function Estimation
In this section, I will be computing L-function via Lest() of spatstat package which is normalises the K-function to a linear scale for easier interpretation.
6.2.1 Yinmarbin District
1) Computing L-function Estimation
L_yinmarbin = Lest(yinmarbin_ppp_owin, correction = "Ripley")
plot(L_yinmarbin, . -r ~ r, ylab= "L(d)-r", xlab = "d(km)",
main = paste("Yinmarbin District (L-Function)"))
Observations
The theoretical L-function under CSR remains at 0, meaning that if the points were randomly distributed, the L-function would be expected to show no significant clustering or dispersion.
There is a consistent deviation of each observed L value from the theoretical value across all distances, particularly with the strongest deviation at <5 km and > 20km. This confirms that conflict points are both localised in towns of Yinmarbin and found across large areas in Yinmarbin.
2) Performing Complete Spatial Randomness Test
To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:
- Ho = The distribution of conflict events in Myanmar are randomly distributed.
- H1= The distribution of conflict events in Myanmar are not randomly distributed.
- The null hypothesis will be rejected if the observed L-function lies above/below the theoretical L-function and envelope.
I will also perform monta carlo simulation test using envelope() of the spatstat package.
# Monte Carlo test with L-function
L_yinmarbin_csr <- envelope(yinmarbin_ppp_owin, Lest,
nsim = 39, rank = 1, glocal=TRUE)
plot(L_yinmarbin_csr, . - r ~ r, xlab="d(km)", ylab="L(d)-r",
main = paste("Yinmarbin District (CSR)"))
Observations
The envelope outputted is narrow which once again reflects lower variability in spatial distribution in Yinmarbin where a large number of points found in this district could have led to more stable results. There is also statistically significant clustering particularly in specific towns in Yinmarbin and across Yinmarbin itself, since we see larger deviations of observed L values from the envelope.
6.2.2 Shwebo District
1) Computing L-function Estimation
L_shwebo = Lest(shwebo_ppp_owin, correction = "Ripley")
plot(L_shwebo, . -r ~ r, ylab= "L(d)-r", xlab = "d(km)",
main = paste("Shwebo District (L-Function)"))
Observations
The observed L values obtained for Shwebo districct tend to show a smaller deviation above the theoretical L-function in smaller areas like towns, suggesting localised clustering is less pronounced. Conversely, larger deviations are seen at larger distances which shows broader patterns of clustering are present in Shwebo district.
2) Performing Complete Spatial Randomness Test
# Monte Carlo test with L-function
L_shwebo_csr <- envelope(shwebo_ppp_owin, Lest,
nsim = 39, rank = 1, glocal=TRUE)
plot(L_shwebo_csr, . - r ~ r, xlab="d(km)", ylab="L(d)-r",
main = paste("Shwebo District (CSR)"))
Observations
The envelope size continues to grow larger as distance increases, which indicates greater variability in spatial patterns of the Shwebo district but overall, we can confirm strong clustering of conflicts especially in bigger areas of the district which suggests wide-scale hotspots, as well as some localised clustering at smaller scales.
6.2.3 Pakokku District
1) Computing L-function Estimation
L_pakokku = Lest(pakokku_ppp_owin, correction = "Ripley")
plot(L_yinmarbin, . -r ~ r, ylab= "L(d)-r", xlab = "d(km)",
main = paste("Pakokku District (L-Function)"))
Observations
Like the other districts, there is clustering seen in Pakokku across smaller areas like towns and the distribution across the district itself.
2) Performing Complete Spatial Randomness Test
# Monte Carlo test with L-function
L_pakokku_csr <- envelope(pakokku_ppp_owin, Lest,
nsim = 39, rank = 1, glocal=TRUE)
Observations
Under the CSR test, the intensity of clustering weakens across larger distances (>20km) than we would have assumed in the previous step. This happens as Monte Carlo simulation accounts for random fluctuations especially at larger distances (e.g., >20 km) where random variability in the point pattern naturally increases.
In short, significant clustering is present at both small and large scales, but with diminishing intensity as we measure larger distances.
6.2.4 Mandalay District
1) Computing L-function Estimation
L_mandalay = Lest(mandalay_ppp_owin, correction = "Ripley")
plot(L_mandalay, . -r ~ r, ylab= "L(d)-r", xlab = "d(km)",
main = paste("Mandalay District (L-Function)"))
Observations
The distance scale for Mandalay is the smallest (0-10km) since the district itself has a smaller area. Nonetheless, we see almost a consistent spread of localised and widespread distribution of conflict points in Mandalay which indicates that conflicts occur across both specific towns and the across the districct.
2) Performing Complete Spatial Randomness Test
# Monte Carlo test with L-function
L_mandalay_csr <- envelope(mandalay_ppp_owin, Lest,
nsim = 39, rank = 1, glocal=TRUE)
plot(L_mandalay_csr, . - r ~ r, xlab="d(km)", ylab="L(d)-r",
main = paste("Mandalay District (CSR)"))
Observations
The Monte Carlo simulation shows that there’s statistically significant clustering across all distances since the observed L values consistently lie above the envelope. There is also a marginal increase in variability of random points observed which is inherent in larger distance scales as seen from the slight widening of the envelope as distance increases.